Nothing
######################
#Begin CCTpack
######################
.onLoad <- function(libname, pkgname) {
assign("pkg_globals", new.env(), envir=parent.env(environment()))
}
######################
#Initializations
######################
options(warn=-3)
### Some users have reported memory allocation errors when a high limit is not set here.
suppressMessages(try(memory.limit(10000),silent=TRUE))
suppressMessages(try(memory.limit(20000),silent=TRUE))
options(warn=0)
######################
#Supplementary Functions
######################
Mode <- function(x) {ux <- unique(x); ux[which.max(tabulate(match(x, ux)))] }
probit <- function(x) {probval <- qnorm(x,0,1); return(probval)}
invprobit <- function(x) {invprobval <- pnorm(x,0,1); return(invprobval)}
######################
#CCTpack GUI, invoked with the command cctgui()
######################
cctgui <- function(){
######################
#Sets Default GUI Variables and Frames
######################
guidat <- list();
guidat$tt <- tktoplevel(bg="#CDEED6")
guidat$datframe <- tkframe(guidat$tt, borderwidth = 0,bg="#CDEED6")
guidat$datframe2 <- tkframe(guidat$tt, borderwidth = 0,bg="#CDEED6")
guidat$datframe3 <- tkframe(guidat$tt, borderwidth = 0,bg="#CDEED6")
guidat$datframe4 <- tkframe(guidat$tt, borderwidth = 0,bg="#CDEED6")
guidat$applyframe1 <- tkframe(guidat$tt, borderwidth = 0,bg="#CDEED6")
guidat$applyframe2 <- tkframe(guidat$tt, borderwidth = 0,bg="#CDEED6")
guidat$applyframe3 <- tkframe(guidat$tt, borderwidth = 0,bg="#CDEED6")
guidat$resultsframe1 <- tkframe(guidat$tt, borderwidth = 0,bg="#CDEED6")
guidat$resultsframe2 <- tkframe(guidat$tt, borderwidth = 0,bg="#CDEED6")
guidat$resultsframe3 <- tkframe(guidat$tt, borderwidth = 0,bg="#CDEED6")
guidat$settingsframe <- tkframe(guidat$tt, borderwidth = 0,bg="#CDEED6")
guidat$polywind <- FALSE
guidat$mval <- FALSE
guidat$samplesvar <- tclVar("10000")
guidat$chainsvar <- tclVar("3")
guidat$burninvar <- tclVar("2000")
guidat$thinvar <- tclVar("1")
guidat$culturesvar <- tclVar("1")
guidat$paravar <- tclVar("1")
guidat$itemdiffvar <- tclVar("0")
guidat$polyvar <- tclVar("0")
guidat$varnametry <- tclVar("cctfit")
tkwm.title(guidat$tt,"CCT Model Application Software")
######################
#The GUI Grid Setup
######################
guidat$samples.entry <- tkentry(guidat$settingsframe, textvariable=guidat$samplesvar,width="6")
guidat$chains.entry <- tkentry(guidat$settingsframe, textvariable=guidat$chainsvar,width="2")
guidat$burnin.entry <- tkentry(guidat$settingsframe, textvariable=guidat$burninvar,width="6")
guidat$thin.entry <- tkentry(guidat$settingsframe, textvariable=guidat$thinvar,width="2")
guidat$loaddata.but <- tkbutton(guidat$datframe, text="Load Data", command=loadfilefuncbutton)
guidat$screeplot.but <- tkbutton(guidat$datframe, text="Scree Plot", command=screeplotfuncbutton)
guidat$poly.but <- tkcheckbutton(guidat$datframe4,variable=guidat$polyvar,text="Use polychoric correlations",bg="#CDEED6")
guidat$para.but <- tkcheckbutton(guidat$applyframe2,variable=guidat$paravar,text="Parallel run",bg="#CDEED6")
guidat$item.but <- tkcheckbutton(guidat$applyframe2,variable=guidat$itemdiffvar,text="Item difficulty",bg="#CDEED6")
guidat$cultdown.but <- tkbutton(guidat$applyframe2, text="<", command=cultdownfuncbutton)
guidat$cultup.but <- tkbutton(guidat$applyframe2, text=">", command=cultupfuncbutton)
guidat$applymodel.but <- tkbutton(guidat$applyframe2, text="Apply CCT Model", command=applymodelfuncbutton)
guidat$cultures.entry <- tkentry(guidat$applyframe2, textvariable=guidat$culturesvar,width="2",disabledforeground="black",disabledbackground="white")
guidat$varname.entry <- tkentry(guidat$applyframe3, textvariable=guidat$varnametry,width="10")
tkbind(guidat$varname.entry, "<Key>", valid_inputkey)
tkbind(guidat$varname.entry, "<FocusOut>", valid_inputfo)
tkbind(guidat$varname.entry, "<Return>", valid_inputfo)
guidat$plotresults.but <- tkbutton(guidat$resultsframe1, text = "Plot Results", command = plotresultsfuncbutton)
guidat$doppc.but <- tkbutton(guidat$resultsframe1, text = "Run Checks", command = ppcfuncbutton)
guidat$exportresults.but <- tkbutton(guidat$resultsframe1, text = "Export Results", command = exportfuncbutton)
guidat$memb.but <- tkbutton(guidat$resultsframe2, text = "Cluster Membs", command = membfuncbutton)
guidat$mvest.but <- tkbutton(guidat$resultsframe2, text = "NA Value Est", command = mvestfuncbutton)
guidat$traceplotdiscrete.but <- tkbutton(guidat$resultsframe2, text = "Traceplot Discrete", command = traceplotdiscretefuncbutton)
guidat$traceplotall.but <- tkbutton(guidat$resultsframe2, text = "Traceplot All", command = traceplotallfuncbutton)
guidat$printfit.but <- tkbutton(guidat$resultsframe3, text = "Print Fit", command = printfitfuncbutton)
guidat$summary.but <- tkbutton(guidat$resultsframe3, text = "Fit Summary", command = summaryfuncbutton)
guidat$sendconsole.but <- tkbutton(guidat$resultsframe3, text = "Send Console", command = sendconsolefuncbutton)
guidat$datafiletxt <- tktext(guidat$tt,bg="white",width=20,height=1)
guidat$resptxt <- tktext(guidat$tt,bg="white",width=4,height=1)
guidat$itemtxt <- tktext(guidat$tt,bg="white",width=4,height=1)
guidat$dattypetxt <- tktext(guidat$tt,bg="white",width=11,height=1)
guidat$modeltxt <- tktext(guidat$tt,bg="white",width=4,height=1)
tkgrid(guidat$datframe,columnspan=3,row=1,column=0)
tkgrid(guidat$datframe2,columnspan=4,row=2,column=0)
tkgrid(guidat$datframe3,columnspan=4,row=3,column=0)
tkgrid(guidat$datframe4,columnspan=5,row=4,column=0)
tkgrid(guidat$applyframe1,columnspan=10,row=5,column=0)
tkgrid(guidat$applyframe2,columnspan=10,row=6,column=0)
tkgrid(guidat$applyframe3,columnspan=10,row=7,column=0)
tkgrid(guidat$resultsframe1,columnspan=3,row=8)
tkgrid(guidat$resultsframe2,columnspan=4,row=9)
tkgrid(guidat$resultsframe3,columnspan=3,row=10)
tkgrid(guidat$settingsframe,columnspan=8,row=11)
tkgrid(tklabel(guidat$datframe,text="",bg="#CDEED6"),columnspan=3, pady = 0)
tkgrid(tklabel(guidat$datframe,text="Data Input",bg="#CDEED6"),columnspan=3, pady = 2)
tkgrid(guidat$loaddata.but,guidat$datafiletxt,guidat$screeplot.but,pady= 10, padx= 10)
tkgrid(tklabel(guidat$datframe2,text="Number of Respondents",bg="#CDEED6"),guidat$resptxt,tklabel(guidat$datframe2,text="Number of Items",bg="#CDEED6"),guidat$itemtxt, padx = 2, pady = 5)
tkgrid(tklabel(guidat$datframe3,text="Data Type Detected",bg="#CDEED6"),guidat$dattypetxt,tklabel(guidat$datframe3,text="CCT Model",bg="#CDEED6"),guidat$modeltxt, padx = 2, pady = 5)
tkgrid(tklabel(guidat$applyframe1,text="Model Application",bg="#CDEED6"),columnspan=3, pady = 5)
tkgrid(guidat$cultdown.but, guidat$cultures.entry, guidat$cultup.but, tklabel(guidat$applyframe2,text="Cultures",bg="#CDEED6"),guidat$item.but,guidat$para.but,guidat$applymodel.but,pady= 10, padx= 2)
tkgrid(guidat$varname.entry,tklabel(guidat$applyframe3,text="Name of Fit",bg="#CDEED6"),padx= 2)
tkconfigure(guidat$screeplot.but, state="disabled")
tkgrid(tklabel(guidat$resultsframe1,text="Application Results",bg="#CDEED6"),columnspan=3, pady = 5)
tkgrid(guidat$doppc.but,guidat$plotresults.but,guidat$exportresults.but,pady= 5, padx= 10)
tkgrid(guidat$memb.but,guidat$mvest.but,guidat$traceplotdiscrete.but,guidat$traceplotall.but,pady= 5, padx= 5)
tkgrid(guidat$printfit.but,guidat$summary.but,guidat$sendconsole.but,pady= 5, padx= 10)
tkconfigure(guidat$applymodel.but, state="disabled")
tkconfigure(guidat$plotresults.but, state="disabled")
tkconfigure(guidat$doppc.but, state="disabled")
tkconfigure(guidat$exportresults.but, state="disabled")
tkconfigure(guidat$printfit.but, state="disabled")
tkconfigure(guidat$summary.but, state="disabled")
tkconfigure(guidat$sendconsole.but, state="disabled")
tkconfigure(guidat$memb.but, state="disabled")
tkconfigure(guidat$mvest.but, state="disabled")
tkconfigure(guidat$traceplotdiscrete.but, state="disabled")
tkconfigure(guidat$traceplotall.but, state="disabled")
tkgrid(tklabel(guidat$settingsframe,text="Sampler Settings (Optional)",bg="#CDEED6"),columnspan=8, pady = 5)
tkgrid(tklabel(guidat$settingsframe,text="Samples",bg="#CDEED6"), guidat$samples.entry, tklabel(guidat$settingsframe,text="Chains",bg="#CDEED6"), guidat$chains.entry, tklabel(guidat$settingsframe,text="Burn-in",bg="#CDEED6"), guidat$burnin.entry, tklabel(guidat$settingsframe,text="Thinning",bg="#CDEED6"), guidat$thin.entry, pady= 10, padx= 2)
tkgrid(tklabel(guidat$settingsframe,text="",bg="#CDEED6"),columnspan=8, pady = 0)
tkconfigure(guidat$cultures.entry, bg="white", state="disabled", background="black")
tkinsert(guidat$datafiletxt,"end","(load data file)")
tkconfigure(guidat$datafiletxt, state="disabled")
tkconfigure(guidat$resptxt, state="disabled")
tkconfigure(guidat$itemtxt, state="disabled")
tkconfigure(guidat$dattypetxt, state="disabled")
tkconfigure(guidat$modeltxt, state="disabled")
tkconfigure(guidat$varname.entry, state="disabled")
tkconfigure(guidat$para.but, state="disabled")
tkconfigure(guidat$item.but, state="disabled")
tkconfigure(guidat$cultdown.but, state="disabled")
tkconfigure(guidat$cultup.but, state="disabled")
tkgrid.columnconfigure(guidat$tt,0,weight=1)
tkgrid.rowconfigure(guidat$tt,0,weight=1)
tkwm.resizable(guidat$tt,0,0)
guidat$guidatnames <- names(guidat)
list2env(guidat,pkg_globals)
assign("guidat", guidat, pkg_globals)
message("\n ...Starting CCT Inference Software\n")
}
######################
#Standalone Functions to use instead of the GUI
######################
######################
#All-in-one function to do each of the GUI tasks, if the gui is not compatible for the user, or not preferred
#- Takes in a data as a respondent by item array or matrix
#- 'clusters' defines # of clusters, and 'itemdiff' if item difficulty is wanted
#- 'jags' defines for the JAGS MCMC, the number of samples, chains, burn-in, thinning
#- 'runchecks' if one wants the posterior predictive checks calculated after inference
#- 'exportfilename' set a name different from "" if one wants to automatically export the results
# to the working directory
######################
cctapply <- function(data,clusters=1,itemdiff=FALSE,biases=TRUE,samples=10000,chains=3,burnin=2000,thinning=1,runchecks=FALSE,exportfilename="",polych=FALSE,parallel=FALSE,seed=NULL,plotr=FALSE){
if(!is.null(seed)){
set.seed(1); rm(list=".Random.seed", envir=globalenv())
set.seed(seed)
}else{set.seed(1); rm(list=".Random.seed", envir=globalenv())}
datob <- loadfilefunc(data)
datob <- screeplotfunc(datob,noplot=TRUE)
cctfit <- applymodelfunc(datob,clusters=clusters,itemdiff=itemdiff,biases=biases,jags.iter=samples,jags.chains=chains,jags.burnin=burnin,jags.thin=thinning,parallel=parallel)
if(plotr==TRUE){plotresultsfunc(cctfit)}
if(runchecks==TRUE){
cctfit <- ppcfunc(cctfit,polych=polych,doplot=FALSE)
}
if(exportfilename!=""){
exportfunc(cctfit,filename=exportfilename)
}
return(cctfit)
}
######################
#Manual function to do the scree plot, equivalent to the 'Scree Plot Button'
#- Takes in a data as a respondent by item array or matrix
######################
cctscree <- function(data,polych=FALSE){
datob <- loadfilefunc(data)
datob <- screeplotfunc(datob,polych=polych)
}
######################
#Manual function to plot the results, equivalent to the 'Plot Results Button'
#- Takes the cctfit object from cctapply() or the 'Apply CCT Model' button
######################
cctresults <- function(cctfit){
plotresultsfunc(cctfit)
}
######################
#Manual function to plot the posterior predictive check plots, equivalent to the 'Run Checks Button'
#- Takes the cctfit object from cctapply() or the 'Apply CCT Model' button
#- Plots the posterior predictive checks, or calculates them and then plots them (if they have not been calculated yet)
######################
cctchecks <- cctppc <- function(cctfit,polych=FALSE,doplot=TRUE){
if(cctfit$whmodel=="LTRM" && cctfit$checksrun==TRUE){
if(cctfit$polycor!=polych){cctfit$checksrun <- FALSE}
}
cctfit <- ppcfunc(cctfit,polych=polych)
return(cctfit)
}
######################
#Manual function to export the results, equivalent to the 'Export Results Button'
#- Takes the cctfit object from cctapply() or the 'Apply CCT Model' button
#- Exports the cctfit object as a .Rdata file, and plots of the scree plot,
# posterior results, and posterior predictive checks
######################
cctexport <- function(cctfit,filename="CCTpackdata.Rdata"){
cctfit <- exportfunc(cctfit,filename=filename)
}
######################
#End of Standalone Functions
######################
#######################
#Background Functions for the GUI and/or Standalone Functions
######################
######################
#Function for the 'Load Data' Button
#- Detects the type of data
#- Detects the number of respondents and number of items
#- Selects the appropriate model for the data
#- Detects the number of missing data points (NA values)
#- Estimates initial estimates for missing data (if any) for the initial scree plot
#- Reports this information back to the GUI or R console
#- Enables the Apply Model Button on the GUI
#- Disables buttons GUI buttons: 'Run Checks', 'Plot Results', 'Export Results'
# if they are active from a previous run
######################
loadfilefuncbutton <- function(){
loadfilefunc(gui=TRUE)
}
loadfilefunc <- function(data=0,gui=FALSE,polych=FALSE){
if(gui==TRUE){guidat <- get("guidat", pkg_globals)}
datob <- list(); datob$polych <- polych; datob$datfactorsp <- 1; datob$datfactorsc
if(gui==TRUE){
datob$fileName <- file.path(tclvalue(tkgetOpenFile(filetypes = "{{csv Files} {.csv .txt}}")))
if (!nchar(datob$fileName)) {
return()
}else{
datob$dat <- as.matrix(read.csv(datob$fileName,header=FALSE))
if(is.list(datob$dat)){
datob$dat <- matrix(as.numeric(unlist(datob$dat)),dim(datob$dat)[1],dim(datob$dat)[2])
}
options(warn=-3)
if(!is.numeric(datob$dat[1,1])){datob$dat <- matrix(as.numeric(datob$dat),dim(datob$dat)[1],dim(datob$dat)[2])
if(all(is.na(datob$dat[,1]))){datob$dat <- datob$dat[,-1];
}
if(all(is.na(datob$dat[1,]))){datob$dat <- datob$dat[-1,];
}
}
if(sum(is.na(datob$dat))>0){
datob$thena <- which(is.na(datob$dat),arr.ind=TRUE)
datob$thenalist <- which(is.na(datob$dat))
datob$dat[datob$thena] <- datob$dat[max(which(!is.na(datob$dat)),arr.ind=TRUE)]
datob$mval <- TRUE
}else{datob$mval <- FALSE}
options(warn=0)
if(all(datob$dat[1:dim(datob$dat)[1],1]==1:dim(datob$dat)[1])){datob$dat <- datob$dat[,-1]; if(datob$mval==TRUE){datob$thena[,2] <- datob$thena[,2] - 1}}
setwd(file.path(dirname(datob$fileName)))
tkconfigure(guidat$screeplot.but, state="normal")
tkconfigure(guidat$applymodel.but, state="normal")
}
tkconfigure(guidat$datafiletxt, state="normal")
tkconfigure(guidat$resptxt, state="normal")
tkconfigure(guidat$itemtxt, state="normal")
tkconfigure(guidat$dattypetxt, state="normal")
tkconfigure(guidat$modeltxt, state="normal")
tkconfigure(guidat$para.but, state="normal")
tkconfigure(guidat$item.but, state="normal")
tkconfigure(guidat$cultdown.but, state="normal")
tkconfigure(guidat$cultup.but, state="normal")
tkdelete(guidat$datafiletxt,"1.0","800.0")
tkdelete(guidat$resptxt,"1.0","800.0")
tkdelete(guidat$itemtxt,"1.0","800.0")
tkdelete(guidat$dattypetxt,"1.0","800.0")
tkdelete(guidat$modeltxt,"1.0","800.0")
tkinsert(guidat$datafiletxt,"end",basename(datob$fileName)); datob$nmiss <- 0
is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol
if(!all(is.wholenumber(datob$dat))){
invlogit <- function(x){x <- 1 / (1 + exp(-x)); return(x)}
logit <- function(x){x <- log(x/(1-x)); return(x)}
tkinsert(guidat$modeltxt,"end","CRM")
datob$whmodel <- "CRM"
if(min(datob$dat) >= 0 && max(datob$dat) <= 1){
datob$datatype <- "Continuous"
message("\n ...Continuous data detected")
tkinsert(guidat$dattypetxt,"end","Continuous in [0,1]")
datob$dat[datob$dat<=0] <- .001; datob$dat[datob$dat>=1] <- .999
datob$dat <- logit(datob$dat)
if(datob$mval==TRUE){
datob$dat[datob$thena] <- NA
datob$thenalist <- which(is.na(datob$dat))
datob$dat[datob$thena] <- colMeans(datob$dat[,datob$thena[,2]],na.rm=TRUE)
}
}else{
datob$datatype <- "Continuous"
tkinsert(guidat$dattypetxt,"end","Continuous")
datob$dat <- invlogit(datob$dat)
datob$dat[datob$dat<=0] <- .001; datob$dat[datob$dat>=1] <- .999
datob$dat <- logit(datob$dat)
if(datob$mval==TRUE){
datob$dat[datob$thena] <- NA
datob$thenalist <- which(is.na(datob$dat))
datob$dat[datob$thena] <- colMeans(datob$dat[,datob$thena[,2]],na.rm=TRUE)
}
datob$datatype <- "Continuous"
message("\n ...Continuous data detected")
}
}else{if(min(datob$dat) >= 0 && max(datob$dat) <= 1){
tkinsert(guidat$modeltxt,"end","GCM")
datob$whmodel <- "GCM"
datob$datatype <- "Binary"
tkinsert(guidat$dattypetxt,"end","Binary")
message("\n ...Binary (dichotomous) data detected")
if(datob$mval==TRUE){
datob$dat[datob$thena] <- NA
datob$thenalist <- which(is.na(datob$dat))
for(i in unique(datob$thena[,2])){
datob$dat[,i][is.na(datob$dat[,i])] <- Mode(datob$dat[,i][which(!is.na(datob$dat[,i]))])
}
}
}else{
tkinsert(guidat$modeltxt,"end","LTRM")
datob$whmodel <- "LTRM"
datob$datatype <- "Ordinal"
tkinsert(guidat$dattypetxt,"end","Ordinal")
message("\n ...Ordinal (categorical) data detected")
if(guidat$polywind==FALSE){
tkgrid(guidat$poly.but,pady= 10, padx= 2)
guidat$polywind <- TRUE
}
if(datob$mval==TRUE){
datob$dat[datob$thena] <- NA
datob$thenalist <- which(is.na(datob$dat))
for(i in unique(datob$thena[,2])){
datob$dat[,i][is.na(datob$dat[,i])] <- Mode(datob$dat[,i][which(!is.na(datob$dat[,i]))])
}
}
}
}
tkinsert(guidat$resptxt,"end",dim(datob$dat)[1])
tkinsert(guidat$itemtxt,"end",dim(datob$dat)[2])
tkconfigure(guidat$datafiletxt, state="disabled")
tkconfigure(guidat$resptxt, state="disabled")
tkconfigure(guidat$itemtxt, state="disabled")
tkconfigure(guidat$dattypetxt, state="disabled")
tkconfigure(guidat$modeltxt, state="disabled")
datob$datind <- cbind(expand.grid(t(row(datob$dat))),expand.grid(t(col(datob$dat))),expand.grid(t(datob$dat)))
if(datob$mval==TRUE){
datob$nmiss <- dim(datob$thena)[1]
datob$datind <- datob$datind[-datob$thenalist,]
datob$datna <- datob$dat
datob$datna[datob$thena] <- NA
message("\n ...Data has ",datob$nmiss," missing values out of ",length(datob$dat))
}
if(datob$whmodel=="LTRM"){
tkconfigure(guidat$poly.but, state="normal")
}else{
if(guidat$polywind==TRUE){
tkconfigure(guidat$poly.but, state="disabled")
} }
tkconfigure(guidat$varname.entry, state="normal")
tkconfigure(guidat$plotresults.but, state="disabled")
tkconfigure(guidat$doppc.but, state="disabled")
tkconfigure(guidat$exportresults.but, state="disabled")
tkconfigure(guidat$printfit.but, state="disabled")
tkconfigure(guidat$summary.but, state="disabled")
tkconfigure(guidat$sendconsole.but, state="disabled")
tkconfigure(guidat$memb.but, state="disabled")
tkconfigure(guidat$mvest.but, state="disabled")
tkconfigure(guidat$traceplotdiscrete.but, state="disabled")
tkconfigure(guidat$traceplotall.but, state="disabled")
datob$nobs <- dim(datob$datind)[1]
datob$datobnames <- names(datob)
#list2env(datob,pkg_globals)
assign("datob", datob, pkg_globals)
guidat$n <- dim(datob$dat)[1]; guidat$m <- dim(datob$dat)[2]; guidat$mval <- datob$mval
assign("guidat", guidat, pkg_globals)
message("\n ...Data loaded")
}
if(gui==FALSE){
if(is.list(data)){
datob$dat <- matrix(as.numeric(unlist(data)),dim(data)[1],dim(data)[2])
}else{
datob$dat <- data
}
options(warn=-3)
if(!is.numeric(datob$dat[1,1])){datob$dat <- matrix(as.numeric(datob$dat),dim(datob$dat)[1],dim(datob$dat)[2])
if(all(is.na(datob$dat[,1]))){datob$dat <- datob$dat[,-1];
}
if(all(is.na(datob$dat[1,]))){datob$dat <- datob$dat[-1,];
}
}
datob$nmiss <- 0
if(sum(is.na(datob$dat))>0){
datob$thena <- which(is.na(datob$dat),arr.ind=TRUE)
datob$thenalist <- which(is.na(datob$dat))
datob$dat[datob$thena] <- unlist(datob$dat)[max(which(!is.na(datob$dat)),arr.ind=TRUE)]
datob$mval <- TRUE
}else{datob$mval <- FALSE}
options(warn=0)
if(all(datob$dat[1:dim(datob$dat)[1],1]==1:dim(datob$dat)[1])){datob$dat <- datob$dat[,-1]; if(datob$mval==TRUE){datob$thena[,2] <- datob$thena[,2] - 1}}
is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol
if(!all(is.wholenumber(datob$dat))){
invlogit <- function(x){x <- 1 / (1 + exp(-x)); return(x)}
logit <- function(x){x <- log(x/(1-x)); return(x)}
datob$whmodel <- "CRM"
if(min(datob$dat) >= 0 && max(datob$dat) <= 1){
datob$datatype <- "Continuous"
message("\n ...Continuous data detected")
datob$dat[datob$dat<=0] <- .001; datob$dat[datob$dat>=1] <- .999
#datob$dat <- logit(datob$dat)
datob$dat <- logit(datob$dat)
if(datob$mval==TRUE){
datob$nmiss <- dim(datob$thena)[1]
datob$dat[datob$thena] <- NA
datob$thenalist <- which(is.na(datob$dat))
datob$dat[datob$thena] <- colMeans(datob$dat[,datob$thena[,2]],na.rm=TRUE)
}
}else{
datob$dat <- invlogit(datob$dat)
datob$dat[datob$dat<=0] <- .001; datob$dat[datob$dat>=1] <- .999
datob$dat <- logit(datob$dat)
if(datob$mval==TRUE){
datob$nmiss <- dim(datob$thena)[1]
datob$dat[datob$thena] <- NA
datob$thenalist <- which(is.na(datob$dat))
datob$dat[datob$thena] <- colMeans(datob$dat[,datob$thena[,2]],na.rm=TRUE)
}
message("\n ...Continuous data detected")}
}else{if(min(datob$dat) >= 0 && max(datob$dat) <= 1){
datob$whmodel <- "GCM"
datob$datatype <- "Binary"
message("\n ...Binary (dichotomous) data detected")
if(datob$mval==TRUE){
datob$nmiss <- dim(datob$thena)[1]
datob$dat[datob$thena] <- NA
datob$thenalist <- which(is.na(datob$dat))
for(i in unique(datob$thena[,2])){
datob$dat[,i][is.na(datob$dat[,i])] <- Mode(datob$dat[,i][which(!is.na(datob$dat[,i]))])
}
}
}else{
datob$whmodel <- "LTRM"
datob$datatype <- "Ordinal"
message("\n ...Ordinal (categorical) data detected")
if(datob$mval==TRUE){
datob$nmiss <- dim(datob$thena)[1]
datob$dat[datob$thena] <- NA
datob$thenalist <- which(is.na(datob$dat))
for(i in unique(datob$thena[,2])){
datob$dat[,i][is.na(datob$dat[,i])] <- Mode(datob$dat[,i][which(!is.na(datob$dat[,i]))])
}
}
}
}
message("\n ...",dim(datob$dat)[1]," respondents and ",dim(datob$dat)[2]," items")
datob$datind <- cbind(expand.grid(t(row(datob$dat))),expand.grid(t(col(datob$dat))),expand.grid(t(datob$dat)))
if(datob$mval==TRUE){message("\n ...Data has ",dim(datob$thena)[1]," missing values out of ",length(datob$dat))
datob$datind <- datob$datind[-datob$thenalist,]
datob$datna <- datob$dat
datob$datna[datob$thena] <- NA
}
datob$nobs <- dim(datob$datind)[1]
return(datob)
}
}
######################
#Function for the 'Scree Plot' Button
#- Uses the data object from loaddatafunc
#- Performs factor analysis on the Pearson correlations of the data
# using the fa() function from psych package
#- If polychoric correlations are picked for the LTRM, uses fa() on the polychoric
# correlations instead using polychoric() from the polycor package
#- Creates a plot of 8 eigenvalues with appropriate titles and labels
#- When exporting the results, this function is run and produces .eps and .jpeg's of the plot
######################
screeplotfuncbutton <- function() {
datob <- get("datob", pkg_globals)
guidat <- get("guidat", pkg_globals)
screeplotfunc(datob=datob,saveplots=0,savedir="",polych=as.logical(as.numeric(tclvalue(guidat$polyvar))))
}
screeplotfunc <- function(datob,saveplots=0,savedir="",gui=FALSE,polych=FALSE,noplot=FALSE) {
options(warn=-3)
if(datob$whmodel!="LTRM"){
if(saveplots==0 && noplot==FALSE){
tmp <- ""
if(datob$mval==TRUE && datob$whmodel=="GCM"){tmp <- ", missing data handled by mode of respective columns"}
if(datob$mval==TRUE && datob$whmodel=="CRM"){tmp <- ", missing data handled by mean of respective columns"}
message("\n ...Producing Scree Plot",tmp)
}
if(saveplots==1){jpeg(file.path(gsub(".Rdata","scree.jpg",savedir)),width = 6, height = 6, units = "in", pointsize = 12,quality=100,res=400)}
if(saveplots==2){postscript(file=file.path(gsub(".Rdata","scree.eps",savedir)), onefile=FALSE, horizontal=FALSE, width = 6, height = 6, paper="special", family="Times")}
if(noplot==FALSE){
if(sum(apply(datob$dat,1,function(x) sd(x,na.rm=TRUE))==0) > 0){
tmp <- datob$dat
if(sum(apply(datob$dat,1,function(x) sd(x,na.rm=TRUE))==0)==1){
tmp[which(apply(tmp,1,function(x) sd(x,na.rm=TRUE))==0),1] <- min(tmp[which(apply(tmp,1,function(x) sd(x,na.rm=TRUE))==0),1]+.01,.99)
}else{
tmp[which(apply(tmp,1,function(x) sd(x,na.rm=TRUE))==0),1] <- sapply(apply(tmp[which(apply(tmp,1,function(x) sd(x,na.rm=TRUE))==0),],1,mean),function(x) min(x+.01,.99))
}
par(oma=c(0,0,0,0),mar=c(4,4,3,1),mgp=c(2.25,.75,0),mfrow=c(1,1))
datob$factors <- suppressMessages(fa(cor(t(datob$dat)))$values[1:8])
suppressMessages(plot(datob$factors,las=1,type="b",bg="black",pch=21,xlab="Factor",ylab="Magnitude",main="Scree Plot of Data"))
}else{
par(oma=c(0,0,0,0),mar=c(4,4,3,1),mgp=c(2.25,.75,0),mfrow=c(1,1))
datob$factors <- suppressMessages(fa(cor(t(datob$dat)))$values[1:8])
suppressMessages(plot(datob$factors,las=1,type="b",bg="black",pch=21,xlab="Factor",ylab="Magnitude",main="Scree Plot of Data"))
}
}
if(saveplots==1 || saveplots ==2){dev.off()}
}else{
if(saveplots==0 && noplot==FALSE){
tmp <- ""
if(datob$mval==TRUE && datob$whmodel=="LTRM"){tmp <- ", missing data handled by mode of respective columns"}
message("\n ...Producing Scree Plot",tmp)
}
if(saveplots==1){jpeg(file.path(gsub(".Rdata","scree.jpg",savedir)),width = 6, height = 6, units = "in", pointsize = 12,quality=100,res=400)}
if(saveplots==2){postscript(file=file.path(gsub(".Rdata","scree.eps",savedir)), onefile=FALSE, horizontal=FALSE, width = 6, height = 6, paper="special", family="Times")}
if(polych==TRUE){
if(length(datob$datfactorsp)==1){
message(" note: utilizing polychoric correlations (time intensive)")
datob$factors <- datob$datfactorsp <- suppressMessages(fa(polychoric(t(datob$dat),global=FALSE)$rho)$values[1:8])
}
datob$datfactors <- datob$datfactorsp
}else{
if(length(datob$datfactorsp)==1){
datob$factors <- datob$datfactorsc <- suppressMessages(fa(cor(t(datob$dat)))$values[1:8])
}
datob$datfactors <- datob$datfactorsc
}
if(noplot==FALSE){
par(oma=c(0,0,0,0),mar=c(4,4,3,1),mgp=c(2.25,.75,0),mfrow=c(1,1))
plot(datob$datfactors,las=1,type="b",bg="black",pch=21,xlab="Factor",ylab="Magnitude",main="Scree Plot of Data")
}
if(gui==TRUE){
datob$datobnames <- names(datob)
#list2env(datob,pkg_globals)
assign("datob", datob, pkg_globals)
}
if(saveplots==1 || saveplots ==2){dev.off()}
}
if(gui==FALSE){return(datob)}
options(warn=0)
}
######################
#Functions for the Increasing/Decreasing Cultures Buttons
######################
#For testing
#itemdifffuncbutton <- function(){print(tclvalue(get("guidat", pkg_globals)$itemdiffvar)) }
cultdownfuncbutton <- function(){guidat <- get("guidat", pkg_globals); guidat$culturesvar <- tclVar(as.character(max(1,as.numeric(tclvalue(guidat$culturesvar))-1))); tkconfigure(guidat$cultures.entry, textvariable=guidat$culturesvar); assign("guidat", guidat, pkg_globals)}
cultupfuncbutton <- function(){guidat <- get("guidat", pkg_globals); guidat$culturesvar <- tclVar(as.character(min(floor(guidat$n/4),as.numeric(tclvalue(guidat$culturesvar))+1))); tkconfigure(guidat$cultures.entry, textvariable=guidat$culturesvar); assign("guidat", guidat, pkg_globals)}
#cultdownfuncbutton <- function(){guidat <- get("guidat", pkg_globals); guidat$culturesvar <- tclVar(as.character(max(1,as.numeric(tclvalue(guidat$culturesvar))-1))); tkconfigure(guidat$cultures.entry, textvariable=guidat$culturesvar); assign("guidat", guidat, pkg_globals)}
#cultupfuncbutton <- function(){guidat <- get("guidat", pkg_globals); guidat$culturesvar <- tclVar(as.character(min(floor(dim(datob$dat)[1]/4),as.numeric(tclvalue(guidat$culturesvar))+1))); tkconfigure(guidat$cultures.entry, textvariable=guidat$culturesvar); assign("guidat", guidat, pkg_globals)}
######################
#Functions for Inputting Object Name Entry Field
######################
valid_inputkey <- function() {
guidat <- get("guidat", pkg_globals);
val <- gsub(" ","",tclvalue(guidat$varnametry))
nchars <- nchar(make.names(val))
if(val != ""){
if(nchars > 10){
guidat$varnametry <- tclVar(as.character(substr(make.names(val),start=1,stop=10)))
tkconfigure(guidat$varname.entry, textvariable=guidat$varnametry)
}else{
guidat$varnametry <- tclVar(as.character(make.names(val)))
tkconfigure(guidat$varname.entry, textvariable=guidat$varnametry)
}}
assign("guidat", guidat, pkg_globals)
}
valid_inputfo <- function() {
guidat <- get("guidat", pkg_globals);
val <- gsub(" ","",tclvalue(guidat$varnametry))
nchars <- nchar(make.names(val))
if(val == ""){
guidat$varnametry <- tclVar("cctfit")
tkconfigure(guidat$varname.entry, textvariable=guidat$varnametry)
assign("guidat", guidat, pkg_globals)
return()
}
if(nchars>10){
guidat$varnametry <- tclVar(as.character(substr(make.names(val),start=1,stop=10)))
tkconfigure(guidat$varname.entry, textvariable=guidat$varnametry)
}else{
guidat$varnametry <- tclVar(as.character(make.names(val)))
tkconfigure(guidat$varname.entry, textvariable=guidat$varnametry)
}
assign("guidat", guidat, pkg_globals)
}
######################
#Function for the 'Apply CCT Model' Button
#- Uses the data object from loaddatafunc
#- Applies hierarchical Bayesian inference for the model using JAGS by packages rjags and R2jags
#- Model code for each of the three models are at the bottom of this rcode file
#- Reads in user preferences for model specifications from the GUI
#- Uses these specifications to apply different form(s) of the model
#- During mixture model cases, applies an algorithm that corrects for label-switching/mixing issues
#- Recalculates the statistics for all parameters after correcting for label-switching, as well as
# the Rhat statistics, and DIC
#- Enables the user to look at traceplots for discrete nodes via dtraceplots()
#- Provides the model-based clustering by cctfit$respmem (respondent membership)
#- Provides cctfit$Lik, which is the likelihood of the model evaluated at each sample
#- Enables GUI buttons: 'Run Checks', 'Plot Results', 'Export Results'
######################
applymodelfuncbutton <- function() {
guidat <- get("guidat", pkg_globals)
datob <- get("datob", pkg_globals)
guidat$varname <- tclvalue(guidat$varnametry)
tkconfigure(guidat$varname.entry, state="disabled")
assign("guidat", guidat, pkg_globals)
assign(guidat$varname,
applymodelfunc(datob=datob,clusters=as.numeric(tclvalue(guidat$culturesvar)),itemdiff=as.logical(as.numeric(tclvalue(guidat$itemdiffvar))),
jags.iter=as.numeric(tclvalue(guidat$samplesvar)),jags.chains= as.numeric(tclvalue(guidat$chainsvar)),jags.burnin=as.numeric(tclvalue(guidat$burninvar)),
jags.thin=as.numeric(tclvalue(guidat$thinvar)),parallel= as.logical(as.numeric(tclvalue(guidat$paravar))),gui=TRUE),
pkg_globals
)
#User-friendly (for novices) but not CRAN compliant because it assigns the fit object in the Global environment
#assign(guidat$varname,
#applymodelfunc(datob=datob,clusters=as.numeric(tclvalue(guidat$culturesvar)),itemdiff=as.logical(as.numeric(tclvalue(guidat$itemdiffvar))),
# jags.iter=as.numeric(tclvalue(guidat$samplesvar)),jags.chains= as.numeric(tclvalue(guidat$chainsvar)),jags.burnin=as.numeric(tclvalue(guidat$burninvar)),
# jags.thin=as.numeric(tclvalue(guidat$thinvar)),parallel= as.logical(as.numeric(tclvalue(guidat$paravar))),gui=TRUE),
#inherits=TRUE
#)
}
applymodelfunc <- function(datob,clusters=1,itemdiff=FALSE,biases=TRUE,jags.iter=10000,jags.chains=3,jags.burnin=2000,jags.thin=1,parallel=FALSE,gui=FALSE) {
if(gui==TRUE){guidat <- get("guidat", pkg_globals)}
######################
#Model Codes
#- Used by the 'Apply CCT Model' button/function
#- Each model has 2 code versions, 1 without item difficulty, 1 with item difficulty
######################
mcgcm <-
"model{
for (l in 1:nobs){
D[Y[l,1],Y[l,2]] <- (th[Y[l,1]]*(1-lam[Y[l,2]])) / ((th[Y[l,1]]*(1-lam[Y[l,2]]))+(lam[Y[l,2]]*(1-th[Y[l,1]])))
pY[Y[l,1],Y[l,2]] <- (D[Y[l,1],Y[l,2]]*Z[Y[l,2],Om[Y[l,1]]]) +((1-D[Y[l,1],Y[l,2]])*g[Y[l,1]])
Y[l,3] ~ dbern(pY[Y[l,1],Y[l,2]]) }
for (i in 1:nresp){
Om[i] ~ dcat(pi)
th[i] ~ dbeta(thmu[Om[i]]*thtau[Om[i]],(1-thmu[Om[i]])*thtau[Om[i]])
g[i] ~ dbeta(gmu[Om[i]]*gtau[Om[i]],(1-gmu[Om[i]])*gtau[Om[i]]) }
for (k in 1:nitem){
lam[k] <- .5
for (v in 1:V){
Z[k,v] ~ dbern(p[v]) }}
#Hyper Parameters
gsmu <- 10
gssig <- 10
dsmu <- 10
dssig <- 10
pi[1:V] ~ ddirch(L)
alpha <- 2
for (v in 1:V){
p[v] ~ dunif(0,1)
gmu[v] <- .5
gtau[v] ~ dgamma(pow(gsmu,2)/pow(gssig,2),gsmu/pow(gssig,2))
thmu[v] ~ dbeta(alpha,alpha)
thtau[v] ~ dgamma(pow(dsmu,2)/pow(dssig,2),dsmu/pow(dssig,2))
L[v] <- 1 }
}"
mcgcmid <-
"model{
for (l in 1:nobs){
D[Y[l,1],Y[l,2]] <- (th[Y[l,1]]*(1-lam[Y[l,2],Om[Y[l,1]]])) / ((th[Y[l,1]]*(1-lam[Y[l,2],Om[Y[l,1]]]))+(lam[Y[l,2],Om[Y[l,1]]]*(1-th[Y[l,1]])))
pY[Y[l,1],Y[l,2]] <- (D[Y[l,1],Y[l,2]]*Z[Y[l,2],Om[Y[l,1]]]) +((1-D[Y[l,1],Y[l,2]])*g[Y[l,1]])
Y[l,3] ~ dbern(pY[Y[l,1],Y[l,2]]) }
for (i in 1:nresp){
Om[i] ~ dcat(pi)
th[i] ~ dbeta(thmu[Om[i]]*thtau[Om[i]],(1-thmu[Om[i]])*thtau[Om[i]])
g[i] ~ dbeta(gmu[Om[i]]*gtau[Om[i]],(1-gmu[Om[i]])*gtau[Om[i]]) }
for (k in 1:nitem){
for (v in 1:V){
Z[k,v] ~ dbern(p[v])
lam[k,v] ~ dbeta(lammu[v]*lamtau[v],(1-lammu[v])*lamtau[v])
}}
#Hyper Parameters
lamsmu <- 10
lamssig <- 10
gsmu <- 10
gssig <- 10
dsmu <- 10
dssig <- 10
alpha <- 2
pi[1:V] ~ ddirch(L)
for (v in 1:V){
p[v] ~ dunif(0,1)
lammu[v] <- .5
lamtau[v] ~ dgamma(pow(lamsmu,2)/pow(lamssig,2),lamsmu/pow(lamssig,2))
gmu[v] <- .5
gtau[v] ~ dgamma(pow(gsmu,2)/pow(gssig,2),gsmu/pow(gssig,2))
thmu[v] ~ dbeta(alpha,alpha)
thtau[v] ~ dgamma(pow(dsmu,2)/pow(dssig,2),dsmu/pow(dssig,2))
L[v] <- 1
}
}"
mcgcmidsingle <-
"model{
for (l in 1:nobs){
D[Y[l,1],Y[l,2]] <- (th[Y[l,1]]*(1-lam[Y[l,2],Om[Y[l,1]]])) / ((th[Y[l,1]]*(1-lam[Y[l,2],Om[Y[l,1]]]))+(lam[Y[l,2],Om[Y[l,1]]]*(1-th[Y[l,1]])))
pY[Y[l,1],Y[l,2]] <- (D[Y[l,1],Y[l,2]]*Z[Y[l,2],Om[Y[l,1]]]) +((1-D[Y[l,1],Y[l,2]])*g[Y[l,1]])
Y[l,3] ~ dbern(pY[Y[l,1],Y[l,2]]) }
for (i in 1:nresp){
Om[i] ~ dcat(pi)
th[i] ~ dbeta(thmu[Om[i]]*thtau[Om[i]],(1-thmu[Om[i]])*thtau[Om[i]])
g[i] ~ dbeta(gmu[Om[i]]*gtau[Om[i]],(1-gmu[Om[i]])*gtau[Om[i]]) }
for (k in 1:nitem){
lam[k] ~ dbeta(lammu*lamtau,(1-lammu)*lamtau)
for (v in 1:V){
Z[k,v] ~ dbern(p[v])
}}
#Hyper Parameters
lamsmu <- 10
lamssig <- 10
gsmu <- 10
gssig <- 10
dsmu <- 10
dssig <- 10
alpha <- 2
pi[1:V] ~ ddirch(L)
lammu <- .5
lamtau ~ dgamma(pow(lamsmu,2)/pow(lamssig,2),lamsmu/pow(lamssig,2))
for (v in 1:V){
p[v] ~ dunif(0,1)
gmu[v] <- .5
gtau[v] ~ dgamma(pow(gsmu,2)/pow(gssig,2),gsmu/pow(gssig,2))
thmu[v] ~ dbeta(alpha,alpha)
thtau[v] ~ dgamma(pow(dsmu,2)/pow(dssig,2),dsmu/pow(dssig,2))
L[v] <- 1
}
}"
mcltrm <-
"model{
for (l in 1:nobs){
tau[Y[l,1],Y[l,2]] <- pow(E[Y[l,1]],-2)
pY[Y[l,1],Y[l,2],1] <- pnorm((a[Y[l,1]]*gam[1,Om[Y[l,1]]]) + b[Y[l,1]],T[Y[l,2],Om[Y[l,1]]],tau[Y[l,1],Y[l,2]])
for (c in 2:(C-1)){pY[Y[l,1],Y[l,2],c] <- pnorm((a[Y[l,1]]*gam[c,Om[Y[l,1]]]) + b[Y[l,1]],T[Y[l,2],Om[Y[l,1]]],tau[Y[l,1],Y[l,2]]) - sum(pY[Y[l,1],Y[l,2],1:(c-1)])}
pY[Y[l,1],Y[l,2],C] <- (1 - sum(pY[Y[l,1],Y[l,2],1:(C-1)]))
Y[l,3] ~ dcat(pY[Y[l,1],Y[l,2],1:C])
}
#Parameters
for (i in 1:nresp){
Om[i] ~ dcat(pi)
Elog[i] ~ dnorm(Emu[Om[i]],Etau[Om[i]])
E[i] <- exp(Elog[i])
alog[i] ~ dnorm(amu[Om[i]],atau[Om[i]])T(-2.3,2.3)
a[i] <- exp(alog[i])
b[i] ~ dnorm(bmu[Om[i]],btau[Om[i]]) }
for (k in 1:nitem){
for (v in 1:V){
T[k,v] ~ dnorm(Tmu[v],Ttau[v]) }}
for (v in 1:V){
gam[1:(C-1),v] <- sort(tgam2[1:(C-1),v])
for (c in 1:(C-2)){tgam[c,v] ~ dnorm(0,.1)}
tgam2[1:(C-2),v] <- tgam[1:(C-2),v]
tgam2[C-1,v] <- -sum(tgam[1:(C-2),v]) }
pi[1:V] ~ ddirch(L)
#Hyperparameters
for (v in 1:V){
L[v] <- 1
Tmu[v] ~ dnorm(0,.25)
Tsig[v] ~ dunif(.25,3)
Ttau[v] <- pow(Tsig[v],-2)
Emu[v] ~ dnorm(0,.01)
Etau[v] ~ dgamma(.01, .01)
amu[v] <- 0
atau[v] ~ dgamma(.01, .01)T(.01,)
bmu[v] <- 0
btau[v] ~ dgamma(.01, .01)
}
}"
mcltrmid <-
"model{
for (l in 1:nobs){
tau[Y[l,1],Y[l,2]] <- pow(E[Y[l,1]]*exp(lam[Y[l,2],Om[Y[l,1]]]),-2)
pY[Y[l,1],Y[l,2],1] <- pnorm((a[Y[l,1]]*gam[1,Om[Y[l,1]]]) + b[Y[l,1]],T[Y[l,2],Om[Y[l,1]]],tau[Y[l,1],Y[l,2]])
for (c in 2:(C-1)){pY[Y[l,1],Y[l,2],c] <- pnorm((a[Y[l,1]]*gam[c,Om[Y[l,1]]]) + b[Y[l,1]],T[Y[l,2],Om[Y[l,1]]],tau[Y[l,1],Y[l,2]]) - sum(pY[Y[l,1],Y[l,2],1:(c-1)])}
pY[Y[l,1],Y[l,2],C] <- (1 - sum(pY[Y[l,1],Y[l,2],1:(C-1)]))
Y[l,3] ~ dcat(pY[Y[l,1],Y[l,2],1:C])
}
#Parameters
for (i in 1:nresp){
Om[i] ~ dcat(pi)
Elog[i] ~ dnorm(Emu[Om[i]],Etau[Om[i]])
E[i] <- exp(Elog[i])
alog[i] ~ dnorm(amu[Om[i]],atau[Om[i]])T(-2.3,2.3)
a[i] <- exp(alog[i])
b[i] ~ dnorm(bmu[Om[i]],btau[Om[i]])
}
for (k in 1:nitem){
for (v in 1:V){
T[k,v] ~ dnorm(Tmu[v],Ttau[v])
lam[k,v] ~ dnorm(lammu[v],lamtau[v])T(-2.3,2.3)
} }
for (v in 1:V){
gam[1:(C-1),v] <- sort(tgam2[1:(C-1),v])
for (c in 1:(C-2)){tgam[c,v] ~ dnorm(0,.1)}
tgam2[1:(C-2),v] <- tgam[1:(C-2),v]
tgam2[C-1,v] <- -sum(tgam[1:(C-2),v])
}
pi[1:V] ~ ddirch(L)
#Hyperparameters
for (v in 1:V){
L[v] <- 1
Tmu[v] ~ dnorm(0,.25)
Tsig[v] ~ dunif(.25,3)
Ttau[v] <- pow(Tsig[v],-2)
Emu[v] ~ dnorm(0,.01)
Etau[v] ~ dgamma(.01, .01)
amu[v] <- 0
atau[v] ~ dgamma(.01, .01)T(.01,)
bmu[v] <- 0
btau[v] ~ dgamma(.01, .01)
lammu[v] <- 0
lamsig[v] ~ dunif(.25, 2)
lamtau[v] <- pow(lamsig[v],-2)
}
}"
mccrm <-
"model{
for (l in 1:nobs){
Y[l,3] ~ dnorm((a[Y[l,1]]*T[Y[l,2],Om[Y[l,1]]])+b[Y[l,1]],pow(a[Y[l,1]]*E[Y[l,1]],-2))
}
#Parameters
for (i in 1:nresp){
Om[i] ~ dcat(pi)
Elog[i] ~ dnorm(Emu[Om[i]],Etau[Om[i]])
E[i] <- exp(Elog[i])
alog[i] ~ dnorm(amu[Om[i]],atau[Om[i]])T(-2.3,2.3)
a[i] <- exp(alog[i])
b[i] ~ dnorm(bmu[Om[i]],btau[Om[i]]) }
for (k in 1:nitem){
for (v in 1:V){
T[k,v] ~ dnorm(Tmu[v],Ttau[v])
}}
pi[1:V] ~ ddirch(L)
#Hyperparameters
for (v in 1:V){
L[v] <- 1
Tmu[v] ~ dnorm(0,0.25)
Tsig[v] ~ dunif(.25,3)
Ttau[v] <- pow(Tsig[v],-2)
Emu[v] ~ dnorm(0,.01)
Etau[v] ~ dgamma(.01, .01)
amu[v] <- 0
atau[v] ~ dgamma(.01, .01)T(.01,)
bmu[v] <- 0
btau[v] ~ dgamma(.01, .01)
}
}"
mccrmid <-
"model{
for (l in 1:nobs){
Y[l,3] ~ dnorm((a[Y[l,1]]*T[Y[l,2],Om[Y[l,1]]])+b[Y[l,1]],pow(a[Y[l,1]]*E[Y[l,1]]*exp(lam[Y[l,2],Om[Y[l,1]]]),-2))
}
#Parameters
for (i in 1:nresp){
Om[i] ~ dcat(pi)
Elog[i] ~ dnorm(Emu[Om[i]],Etau[Om[i]])
E[i] <- exp(Elog[i])
alog[i] ~ dnorm(amu[Om[i]],atau[Om[i]])T(-2.3,2.3)
a[i] <- exp(alog[i])
b[i] ~ dnorm(bmu[Om[i]],btau[Om[i]]) }
for (k in 1:nitem){
for (v in 1:V){
T[k,v] ~ dnorm(Tmu[v],Ttau[v])
lam[k,v] ~ dnorm(lammu[v],lamtau[v])T(-2.3,2.3)
}}
pi[1:V] ~ ddirch(L)
#Hyperparameters
for (v in 1:V){
L[v] <- 1
Tmu[v] ~ dnorm(0,.25)
Tsig[v] ~ dunif(.25,3)
Ttau[v] <- pow(Tsig[v],-2)
Emu[v] ~ dnorm(0,.01)
Etau[v] ~ dgamma(.01, .01)
amu[v] <- 0
atau[v] ~ dgamma(.01, .01)T(.01,)
bmu[v] <- 0
btau[v] ~ dgamma(.01, .01)
lammu[v] <- 0
lamsig[v] ~ dunif(.25, 2)
lamtau[v] <- pow(lamsig[v],-2)
}
}"
## Diffuse Settings
# L[v] <- 1
# Tmu[v] ~ dnorm(0,.25)
# Ttau[v] ~ dgamma(.01, .01)
# Emu[v] ~ dnorm(0,.01)
# Etau[v] ~ dgamma(.01, .01)
# amu[v] <- 0
# atau[v] ~ dgamma(.01, .01)T(.01,)
# bmu[v] <- 0
# btau[v] ~ dgamma(.01, .01)
# lammu[v] <- 0
# lamtau[v] ~ dgamma(.01, .01)T(.01,)
######################
#Sets up Parameters, Variables, and Data for JAGS for the model selected
######################
if(datob$whmodel=="GCM"){
Y <- datob$datind; nresp <- dim(datob$dat)[1]; nitem <- dim(datob$dat)[2]; V <- clusters; nobs <- dim(datob$datind)[1]
jags.data <- list("Y","nresp","nitem","V","nobs")
if(itemdiff==FALSE){
model.file <- mcgcm
jags.params <- c("Z","th","g","p","thmu","thtau","gmu","gtau","Om","pi")
if(clusters>1){
if(biases==TRUE){
jags.inits <- function(){ list("Z"=matrix(rbinom(nitem*V,1,.5),nitem,V),"th"= runif(nresp,.2,.8), "g"= runif(nresp,.2,.8),"Om"= sample(1:V,nresp,replace=TRUE) )}
}else{
jags.inits <- function(){ list("Z"=matrix(rbinom(nitem*V,1,.5),nitem,V),"th"= runif(nresp,.2,.8), "Om"= sample(1:V,nresp,replace=TRUE) )}
}
}else{
if(biases==TRUE){
jags.inits <- function(){ list("Z"=matrix(rbinom(nitem*V,1,.5),nitem,V),"th"= runif(nresp,.2,.8), "g"= runif(nresp,.2,.8) )}
}else{
jags.inits <- function(){ list("Z"=matrix(rbinom(nitem*V,1,.5),nitem,V),"th"= runif(nresp,.2,.8) )}
}
}
}
if(itemdiff==TRUE){
model.file <- mcgcmid
jags.params <- c("Z","th","g","lam","p","thmu","thtau","gmu","gtau","lammu","lamtau","Om","pi")
if(clusters>1){
if(biases==TRUE){
jags.inits <- function(){ list("Z"=matrix(rbinom(nitem*V,1,.5),nitem,V),"th"= runif(nresp,.2,.8), "g"= runif(nresp,.2,.8), "lam"= matrix(runif(nitem*V,.2,.8),nitem,V), "Om"= sample(1:V,nresp,replace=TRUE) )}
}else{
jags.inits <- function(){ list("Z"=matrix(rbinom(nitem*V,1,.5),nitem,V),"th"= runif(nresp,.2,.8), "lam"= matrix(runif(nitem*V,.2,.8),nitem,V), "Om"= sample(1:V,nresp,replace=TRUE) )}
}
}else{
if(biases==TRUE){
jags.inits <- function(){ list("Z"=matrix(rbinom(nitem*V,1,.5),nitem,V),"th"= runif(nresp,.2,.8), "g"= runif(nresp,.2,.8), "lam"= matrix(runif(nitem*V,.2,.8),nitem,V) )}
}else{
jags.inits <- function(){ list("Z"=matrix(rbinom(nitem*V,1,.5),nitem,V),"th"= runif(nresp,.2,.8), "lam"= matrix(runif(nitem*V,.2,.8),nitem,V) )}
}
}
}
if(clusters==1){
model.file <- gsub(pattern="pi\\[1\\:V\\] ~ ddirch\\(L\\)", "pi <- 1", model.file)
model.file <- gsub(pattern="Om\\[i\\] ~ dcat\\(pi\\)", "Om\\[i\\] <- 1", model.file)
}
if(biases==FALSE){
model.file <- gsub(pattern="g\\[i\\] ~ dbeta\\(gmu\\[Om\\[i\\]\\]\\*gtau\\[Om\\[i\\]\\],\\(1-gmu\\[Om\\[i\\]\\]\\)\\*gtau\\[Om\\[i\\]\\]\\)", "g\\[i\\] <- .5", model.file)
model.file <- gsub(pattern="gtau\\[v\\] ~ dgamma\\(pow\\(gsmu,2\\)/pow\\(gssig,2\\),gsmu/pow\\(gssig,2\\)\\)", "gtau\\[v\\] <- 0", model.file)
}
}
if(datob$whmodel=="LTRM"){
Y <- datob$datind; nresp <- dim(datob$dat)[1]; nitem <- dim(datob$dat)[2]; V <- clusters; C <- max(datob$datind[,3]); nobs <- dim(datob$datind)[1]
jags.data <- list("Y","nresp","nitem","C","V","nobs")
if(itemdiff==FALSE){
model.file <- mcltrm
jags.params <- c("T","gam","E","a","b","Tmu","Ttau","Emu","Etau","amu","atau","bmu","btau","Om","pi")
if(clusters>1){
jags.inits <- function(){ list("tgam"=matrix(rnorm((C-2)*V,0,1),(C-2),V),"T"=matrix(rnorm(nitem*V,0,1),nitem,V),"Emu"= runif(V,.8,2),"Esig"= runif(V,.4,.6), "asig"= runif(V,.4,.6),"bsig"= runif(V,.4,.6),"Om"= sample(1:V,nresp,replace=TRUE) )}
}else{
jags.inits <- function(){ list("tgam"=matrix(rnorm((C-2)*V,0,1),(C-2),V),"T"=matrix(rnorm(nitem*V,0,1),nitem,V),"Emu"= runif(V,.8,2),"Esig"= runif(V,.4,.6), "asig"= runif(V,.4,.6),"bsig"= runif(V,.4,.6) )}
}
if(C==2){
if(clusters>1){
jags.inits <- function(){ list("T"=matrix(rnorm(nitem*V,0,1),nitem,V),"Emu"= runif(V,.8,2),"Esig"= runif(V,.4,.6),"bsig"= runif(V,.4,.6),"Om"= sample(1:V,nresp,replace=TRUE) )}
}else{
jags.inits <- function(){ list("T"=matrix(rnorm(nitem*V,0,1),nitem,V),"Emu"= runif(V,.8,2),"Esig"= runif(V,.4,.6),"bsig"= runif(V,.4,.6) )}
}
model.file <- gsub(pattern="a\\[i\\] <- exp\\(alog\\[i\\]\\)", "a\\[i\\] <- 1", model.file)
model.file <- gsub(pattern="alog\\[i\\] ~ dnorm\\(amu\\[Om\\[i\\]\\],atau\\[Om\\[i\\]\\]\\)T\\(-2.3,2.3\\)", "", model.file)
model.file <- gsub(pattern="atau\\[v\\] ~ dgamma\\(.01,.01\\)T\\(.01,\\)", "atau\\[v\\] <- 0", model.file)
model.file <- gsub(pattern="for \\(c in 2\\:\\(C-1\\)\\)", " #for \\(c in 2\\:\\(C-1\\)\\)", model.file)
model.file <- gsub(pattern="gam\\[1\\:\\(C-1\\),v\\] <- sort\\(tgam2\\[1\\:\\(C-1\\),v\\]\\)", "gam\\[1,v\\] <- 0 }", model.file)
model.file <- gsub(pattern="for \\(c in 1\\:\\(C-2\\)\\)", "#for \\(c in 1\\:\\(C-2\\)\\)", model.file)
model.file <- gsub(pattern="\\(1 - sum\\(pY\\[i,k,1\\:\\(C-1\\)\\]\\)\\)", "1 - pY\\[i,k,1\\]", model.file)
model.file <- gsub(pattern="tgam2\\[1", "#tgam2\\[1", model.file)
model.file <- gsub(pattern="tgam2\\[C", "#tgam2\\[C", model.file)
}
}
if(itemdiff==TRUE){
model.file <- mcltrmid
jags.params <- c("T","lam","gam","E","a","b","Tmu","Ttau","Emu","Etau","amu","atau","bmu","btau","lammu","lamtau","Om","pi")
if(clusters>1){
jags.inits <- function(){ list("tgam"=matrix(rnorm((C-2)*V,0,1),(C-2),V),"T"=matrix(rnorm(nitem*V,0,1),nitem,V),"Emu"= runif(V,.8,2),"Esig"= runif(V,.4,.6), "asig"= runif(V,.4,.6),"bsig"= runif(V,.4,.6),"lamsig"= runif(V,.25,.30), "Om"= sample(1:V,nresp,replace=TRUE) )}
}else{
jags.inits <- function(){ list("tgam"=matrix(rnorm((C-2)*V,0,1),(C-2),V),"T"=matrix(rnorm(nitem*V,0,1),nitem,V),"Emu"= runif(V,.8,2),"Esig"= runif(V,.4,.6), "asig"= runif(V,.4,.6),"bsig"= runif(V,.4,.6),"lamsig"= runif(V,.25,.30) )}
}
if(C==2){
if(clusters>1){
jags.inits <- function(){ list("T"=matrix(rnorm(nitem*V,0,1),nitem,V),"Emu"= runif(V,.8,2),"Esig"= runif(V,.4,.6),"bsig"= runif(V,.4,.6),"lamsig"= runif(V,.25,.30), "Om"= sample(1:V,nresp,replace=TRUE) )}
}else{
jags.inits <- function(){ list("T"=matrix(rnorm(nitem*V,0,1),nitem,V),"Emu"= runif(V,.8,2),"Esig"= runif(V,.4,.6),"bsig"= runif(V,.4,.6),"lamsig"= runif(V,.25,.30) )}
}
model.file <- gsub(pattern="a\\[i\\] <- exp\\(alog\\[i\\]\\)", "a\\[i\\] <- 1", model.file)
model.file <- gsub(pattern="alog\\[i\\] ~ dnorm\\(amu\\[Om\\[i\\]\\],atau\\[Om\\[i\\]\\]\\)T\\(-2.3,2.3\\)", "", model.file)
model.file <- gsub(pattern="atau\\[v\\] ~ dgamma\\(.01, .01\\)T\\(.01,\\)", "atau\\[v\\] <- 0", model.file)
model.file <- gsub(pattern="for \\(c in 2\\:\\(C-1\\)\\)", " #for \\(c in 2\\:\\(C-1\\)\\)", model.file)
model.file <- gsub(pattern="gam\\[1\\:\\(C-1\\),v\\] <- sort\\(tgam2\\[1\\:\\(C-1\\),v\\]\\)", "gam\\[1,v\\] <- 0 }", model.file)
model.file <- gsub(pattern="for \\(c in 1\\:\\(C-2\\)\\)", "#for \\(c in 1\\:\\(C-2\\)\\)", model.file)
model.file <- gsub(pattern="\\(1 - sum\\(pY\\[i,k,1\\:\\(C-1\\)\\]\\)\\)", "1 - pY\\[i,k,1\\]", model.file)
model.file <- gsub(pattern="tgam2\\[1", "#tgam2\\[1", model.file)
model.file <- gsub(pattern="tgam2\\[C", "#tgam2\\[C", model.file)
}
}
if(clusters==1){
model.file <- gsub(pattern="pi\\[1\\:V\\] ~ ddirch\\(L\\)", "pi <- 1", model.file)
model.file <- gsub(pattern="Om\\[i\\] ~ dcat\\(pi\\)", "Om\\[i\\] <- 1", model.file)
}
if(biases==FALSE){
model.file <- gsub(pattern="a\\[i\\] <- exp\\(alog\\[i\\]\\)", "a\\[i\\] <- 1", model.file)
model.file <- gsub(pattern="alog\\[i\\] ~ dnorm\\(amu\\[Om\\[i\\]\\],atau\\[Om\\[i\\]\\]\\)T\\(-2.3,2.3\\)", "alog\\[i\\] <- 0", model.file)
model.file <- gsub(pattern="atau\\[v\\] ~ dgamma\\(.01, .01\\)T\\(.01,\\)", "atau\\[v\\] <- 0", model.file)
model.file <- gsub(pattern="b\\[i\\] ~ dnorm\\(bmu\\[Om\\[i\\]\\],btau\\[Om\\[i\\]\\]) ", "b\\[i\\] <- 0", model.file)
model.file <- gsub(pattern="btau\\[v\\] ~ dgamma\\(.01, .01\\)", "btau\\[v\\] <- 0", model.file)
}
}
if(datob$whmodel=="CRM"){
Y <- datob$datind; nresp <- dim(datob$dat)[1]; nitem <- dim(datob$dat)[2]; V <- clusters; nobs <- dim(datob$datind)[1]
jags.data <- list("Y","nresp","nitem","V","nobs")
if(itemdiff==FALSE){
model.file <- mccrm
jags.params <- c("T","E","a","alog","b","Tmu","Ttau","Emu","Etau","amu","atau","bmu","btau","Om","pi")
if(clusters>1){
jags.inits <- function(){ list("T"=matrix(rnorm(nitem*V,0,.5),nitem,V),"Emu"= runif(V,.8,2),"Esig"= runif(V,.5,.8), "asig"= runif(V,.4,.6),"bsig"= runif(V,.4,.6), "Om"= sample(1:V,nresp,replace=TRUE) )}
}else{
jags.inits <- function(){ list("T"=matrix(rnorm(nitem*V,0,1),nitem,V),"Emu"= runif(V,.8,2),"Esig"= runif(V,.5,.8), "asig"= runif(V,.4,.6),"bsig"= runif(V,.4,.6) )}
}
}
if(itemdiff==TRUE){
model.file <- mccrmid
jags.params <- c("T","E","a","alog","b","lam","Tmu","Ttau","Emu","Etau","amu","atau","bmu","btau","lammu","lamtau","Om","pi")
if(clusters>1){
jags.inits <- function(){ list("T"=matrix(rnorm(nitem*V,0,1),nitem,V),"Emu"= runif(V,.8,2),"Esig"= runif(V,.4,.6), "asig"= runif(V,.4,.6),"bsig"= runif(V,.4,.6),"lamsig"= runif(V,.25,.30), "Om"= sample(1:V,nresp,replace=TRUE) )}
}else{
jags.inits <- function(){ list("T"=matrix(rnorm(nitem*V,0,1),nitem,V),"Emu"= runif(V,.8,2),"Esig"= runif(V,.4,.6), "asig"= runif(V,.4,.6),"bsig"= runif(V,.4,.6),"lamsig"= runif(V,.25,.30) )}
}
}
if(clusters==1){
model.file <- gsub(pattern="pi\\[1\\:V\\] ~ ddirch\\(L\\)", "pi <- 1", model.file)
model.file <- gsub(pattern="Om\\[i\\] ~ dcat\\(pi\\)", "Om\\[i\\] <- 1", model.file)
}
if(biases==FALSE){
model.file <- gsub(pattern="a\\[i\\] <- exp\\(alog\\[i\\]\\)", "a\\[i\\] <- 1", model.file)
model.file <- gsub(pattern="alog\\[i\\] ~ dnorm\\(amu\\[Om\\[i\\]\\],atau\\[Om\\[i\\]\\]\\)T\\(-2.3,2.3\\)", "alog\\[i\\] <- 0", model.file)
model.file <- gsub(pattern="atau\\[v\\] ~ dgamma\\(.01, .01\\)T\\(.01,\\)", "atau\\[v\\] <- 0", model.file)
model.file <- gsub(pattern="b\\[i\\] ~ dnorm\\(bmu\\[Om\\[i\\]\\],btau\\[Om\\[i\\]\\]) ", "b\\[i\\] <- 0", model.file)
model.file <- gsub(pattern="btau\\[v\\] ~ dgamma\\(.01, .01\\)", "btau\\[v\\] <- 0", model.file)
}
}
######################
#Runs the Model in JAGS
#- saves the data object from loadfilefunc within the jags object
#- saves other useful details used later into the jags object
######################
if(parallel==FALSE){
cctfit <- jags(data=jags.data, inits=jags.inits, parameters.to.save=jags.params,
n.chains=jags.chains, n.iter=(jags.iter+jags.burnin), n.burnin=jags.burnin,
n.thin=jags.thin, model.file=textConnection(model.file)) #textConnection(model.file)
cctfit$dataind <- datob$datind; cctfit$data <- datob$dat; cctfit$n <- nresp; cctfit$m <- nitem; cctfit$V <- V
cctfit$mval <- datob$mval; cctfit$itemdiff <- itemdiff; cctfit$checksrun <- FALSE; cctfit$whmodel <- datob$whmodel; cctfit$datob <- datob
if(cctfit$mval==TRUE){cctfit$datamiss <- datob$datna}
if(cctfit$whmodel=="LTRM"){cctfit$C <- C;cctfit$polycor <- FALSE}
}else{
### Prepare / Evaluate parameters for Parallel JAGS run
jags.burnin <- jags.burnin; jags.iter <- jags.iter; jags.chains <- jags.chains; jags.thin <- jags.thin;
tmpfn=tempfile()
tmpcn=file(tmpfn,"w"); cat(model.file,file=tmpcn); close(tmpcn);
message("\n ...Fitting model in parallel \n note: progress bar currently not available with parallel option")
cctfit <- jags.parallel(data=jags.data, inits=jags.inits, parameters.to.save=jags.params,
n.chains=jags.chains, n.iter=(jags.iter+jags.burnin), n.burnin=jags.burnin,
n.thin=jags.thin, model.file=tmpfn, envir=environment(),jags.seed = abs(.Random.seed[3]))
cctfit$dataind <- datob$datind; cctfit$data <- datob$dat; cctfit$n <- nresp; cctfit$m <- nitem; cctfit$V <- V
cctfit$mval <- datob$mval; cctfit$itemdiff <- itemdiff; cctfit$checksrun <- FALSE; cctfit$whmodel <- datob$whmodel; cctfit$datob <- datob
if(cctfit$mval==TRUE){cctfit$datamiss <- datob$datna}
if(cctfit$whmodel=="LTRM"){cctfit$C <- C}
}
######################
#Function used to calculate the Rhats
######################
Rhat1 <- function(mat) {
m <- ncol(mat)
n <- nrow(mat)
b <- apply(mat,2,mean)
B <- sum((b-mean(mat))^2)*n/(m-1)
w <- apply(mat,2,var)
W <- mean(w)
s2hat <- (n-1)/n*W + B/n
Vhat <- s2hat + B/m/n
covWB <- n /m * (cov(w,b^2)-2*mean(b)*cov(w,b))
varV <- (n-1)^2 / n^2 * var(w)/m +
(m+1)^2 / m^2 / n^2 * 2*B^2/(m-1) +
2 * (m-1)*(n-1)/m/n^2 * covWB
df <- 2 * Vhat^2 / varV
R <- sqrt((df+3) * Vhat / (df+1) / W)
return(R)
}
Rhat <- function(arr) {
dm <- dim(arr)
if (length(dm)==2) return(Rhat1(arr))
if (dm[2]==1) return(NULL)
if (dm[3]==1) return(Rhat1(arr[,,1]))
return(apply(arr,3,Rhat1))
}
######################
#Algorithm that corrects for label switching for the GCM
######################
labelswitchalggcm <- function(cctfit,chnind=0){
cctfit2 <- cctfit
nch <- cctfit2$BUGSoutput$n.chains
nsamp <- cctfit2$BUGSoutput$n.keep
if(nch!=1){
ntruths <- cctfit2$V
truths <- array(NA,c(cctfit2$m,nch,ntruths))
inds <- cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Z")]]
inds <- matrix(inds,cctfit2$m,cctfit2$V)
for(v in 1:cctfit$V){truths[,,v] <- t(apply(cctfit$BUGSoutput$sims.array[ ,,inds[,v]],c(2,3),mean))}
V <- cctfit2$V
chstart <- 1
if(length(chnind)==1){
chstart <- 2
chnind <- array(NA,c(V,nch))
chnind[1:V,1] <- 1:V
for(v in 1:V){
for(ch in chstart:nch){
Tind <- c(1:V)[-chnind[,ch][!is.na(chnind[,ch])]]
if(length(Tind)==0){Tind <- c(1:V)}
chnind[v,ch] <- which(max(cor(truths[1:cctfit2$m, 1, v],truths[1:cctfit2$m, ch, Tind]))==cor(truths[1:cctfit2$m, 1, v],truths[1:cctfit2$m, ch, ]))
}}
}
nsamp <- cctfit$BUGSoutput$n.keep
if(cctfit$itemdiff==TRUE){
inds2 <- cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="lam")]]
inds2 <- matrix(inds2,cctfit2$m,cctfit2$V)
inds <- rbind(inds,
inds2,
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="lammu")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="lamtau")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="thmu")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="thtau")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="gmu")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="gtau")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="p")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="pi")]])
}else{
inds <- rbind(inds,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="thmu")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="thtau")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="gmu")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="gtau")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="p")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="pi")]])
}
einds <- cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]]
tmpeinds <- cctfit$BUGSoutput$sims.array[ ,,einds]
for(v in 1:cctfit$V){
for(ch in chstart:nch){
cctfit2$BUGSoutput$sims.array[ ,ch,inds[,v]] <- cctfit$BUGSoutput$sims.array[ ,ch,inds[,chnind[v,ch]]] #this cctfit2 vs. cctfit difference is intentional
if(chstart==2){cctfit2$BUGSoutput$sims.array[ ,ch,einds][cctfit$BUGSoutput$sims.array[ ,ch,einds]==chnind[v,ch]] <- chnind[v,1]}
if(chstart==1){cctfit2$BUGSoutput$sims.array[ ,ch,einds][cctfit$BUGSoutput$sims.array[ ,ch,einds]==v] <- chnind[v,1]}
}}
}
cctfit2$BUGSoutput$sims.list[["Om"]] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]]], c(nsamp*nch,cctfit$n))
cctfit2$BUGSoutput$mean$Om <- apply(cctfit2$BUGSoutput$sims.list[["Om"]],2,mean)
cctfit2$BUGSoutput$sims.matrix <- array(cctfit2$BUGSoutput$sims.array,c(nsamp*nch,dim(cctfit$BUGSoutput$sims.array)[3])) #this cctfit2 vs. cctfit difference is intentional
cctfit2$BUGSoutput$sims.list[["th"]] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="th")]]], c(nsamp*nch,cctfit$n))
cctfit2$BUGSoutput$sims.list[["g"]] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="g")]]], c(nsamp*nch,cctfit$n))
#new code
cctfit2$BUGSoutput$mean[["th"]] <- apply(cctfit2$BUGSoutput$sims.list[["th"]],2,mean)
cctfit2$BUGSoutput$mean[["g"]] <- apply(cctfit2$BUGSoutput$sims.list[["g"]],2,mean)
if(cctfit$itemdiff==TRUE){cctfit2$BUGSoutput$sims.list[["lam"]] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="lam")]]], c(nsamp*nch,cctfit$m))}
cctfit2$BUGSoutput$sims.list[["Z"]] <- array(NA, c(nsamp*nch,cctfit$m,cctfit$V))
if(cctfit$itemdiff==TRUE){
cctfit2$BUGSoutput$sims.list[["lam"]] <- array(NA, c(nsamp*nch,cctfit$m,cctfit$V))
cctfit2$BUGSoutput$sims.list[["lammu"]] <- array(NA, c(nsamp*nch,cctfit$V))
cctfit2$BUGSoutput$sims.list[["lamtau"]] <- array(NA, c(nsamp*nch,cctfit$V))
}
cctfit2$BUGSoutput$sims.list[["thmu"]] <- array(NA, c(nsamp*nch,cctfit$V))
cctfit2$BUGSoutput$sims.list[["thtau"]] <- array(NA, c(nsamp*nch,cctfit$V))
cctfit2$BUGSoutput$sims.list[["gmu"]] <- array(NA, c(nsamp*nch,cctfit$V))
cctfit2$BUGSoutput$sims.list[["gtau"]] <- array(NA, c(nsamp*nch,cctfit$V))
cctfit2$BUGSoutput$sims.list[["p"]] <- array(NA, c(nsamp*nch,cctfit$V))
cctfit2$BUGSoutput$sims.list[["pi"]] <- array(NA, c(nsamp*nch,cctfit$V))
for(v in 1:V){
cctfit2$BUGSoutput$sims.list[["Z"]][,,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Z")]][1:cctfit$m +((v-1)*cctfit$m)]], c(nsamp*nch,cctfit$m))
if(cctfit$itemdiff==TRUE){
cctfit2$BUGSoutput$sims.list[["lam"]][,,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="lam")]][1:cctfit$m +((v-1)*cctfit$m)]], c(nsamp*nch,cctfit$m))
cctfit2$BUGSoutput$sims.list[["lammu"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="lammu")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$sims.list[["lamtau"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="lamtau")]][v]], c(nsamp*nch))
}
cctfit2$BUGSoutput$sims.list[["thmu"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="thmu")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$sims.list[["thtau"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="thtau")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$sims.list[["gmu"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="gmu")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$sims.list[["gtau"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="gtau")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$sims.list[["p"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="p")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$sims.list[["pi"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="pi")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$mean$Z[,v] <- apply(cctfit2$BUGSoutput$sims.list[["Z"]][,,v],2,mean)
if(cctfit$itemdiff==TRUE){
cctfit2$BUGSoutput$mean$lam[,v] <- apply(cctfit2$BUGSoutput$sims.list[["lam"]][,,v],2,mean)
cctfit2$BUGSoutput$mean$lammu[v] <- mean(cctfit2$BUGSoutput$sims.list[["lammu"]][,v])
cctfit2$BUGSoutput$mean$lamtau[v] <- mean(cctfit2$BUGSoutput$sims.list[["lamtau"]][,v])
}
cctfit2$BUGSoutput$mean$thmu[v] <- mean(cctfit2$BUGSoutput$sims.list[["thmu"]][,v])
cctfit2$BUGSoutput$mean$thtau[v] <- mean(cctfit2$BUGSoutput$sims.list[["thtau"]][,v])
cctfit2$BUGSoutput$mean$gmu[v] <- mean(cctfit2$BUGSoutput$sims.list[["gmu"]][,v])
cctfit2$BUGSoutput$mean$gtau[v] <- mean(cctfit2$BUGSoutput$sims.list[["gtau"]][,v])
cctfit2$BUGSoutput$mean$p[v] <- mean(cctfit2$BUGSoutput$sims.list[["p"]][,v])
cctfit2$BUGSoutput$mean$pi[v] <- mean(cctfit2$BUGSoutput$sims.list[["pi"]][,v])
}
if(nch!=1){
cctfit2$BUGSoutput$summary[,1] <- apply(apply(cctfit2$BUGSoutput$sims.array,c(2,3),mean),2,mean)
cctfit2$BUGSoutput$summary[,2] <- apply(apply(cctfit2$BUGSoutput$sims.array,c(2,3),sd),2,sd)
cctfit2$BUGSoutput$summary[,3:7] <- t(apply(cctfit2$BUGSoutput$sims.matrix,2,function(x) quantile(x,probs=c(.025,.25,.50,.75,.975))))
cctfit2$BUGSoutput$summary[,8] <- Rhat(cctfit2$BUGSoutput$sims.array)
cctfit2$BUGSoutput$summary[,8][is.nan(cctfit2$BUGSoutput$summary[,8])] <- 1.000000
dimnames(cctfit2$BUGSoutput$sims.array) <- list(NULL,NULL,rownames(cctfit$BUGSoutput$summary))
dimnames(cctfit2$BUGSoutput$sims.matrix) <- list(NULL,rownames(cctfit$BUGSoutput$summary))
}else{
cctfit2$BUGSoutput$summary[,1] <- apply(cctfit2$BUGSoutput$sims.array,2,mean)
cctfit2$BUGSoutput$summary[,2] <- apply(cctfit2$BUGSoutput$sims.array,2,sd)
cctfit2$BUGSoutput$summary[,3:7] <- t(apply(cctfit2$BUGSoutput$sims.matrix,2,function(x) quantile(x,probs=c(.025,.25,.50,.75,.975))))
cctfit2$BUGSoutput$summary <- cctfit2$BUGSoutput$summary[,-c(8,9)]
dimnames(cctfit2$BUGSoutput$sims.array) <- list(NULL,NULL,rownames(cctfit$BUGSoutput$summary))
dimnames(cctfit2$BUGSoutput$sims.matrix) <- list(NULL,rownames(cctfit$BUGSoutput$summary))
}
cctfit <- cctfit2; rm(cctfit2)
return(cctfit)
}
######################
#Algorithm that corrects for label switching for the LTRM
######################
labelswitchalgltrm <- function(cctfit,chnind=0){
cctfit2 <- cctfit
nch <- cctfit2$BUGSoutput$n.chains
nsamp <- cctfit2$BUGSoutput$n.keep
if(nch!=1){
ntruths <- cctfit2$V
truths <- array(NA,c(cctfit2$m,nch,ntruths))
inds <- cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="T")]]
inds <- matrix(inds,cctfit2$m,cctfit2$V)
for(v in 1:cctfit$V){truths[,,v] <- t(apply(cctfit$BUGSoutput$sims.array[ ,,inds[,v]],c(2,3),mean))}
V <- cctfit2$V
chstart <- 1
if(length(chnind)==1){
chstart <- 2
chnind <- array(NA,c(V,nch))
chnind[1:V,1] <- 1:V
for(v in 1:V){
for(ch in chstart:nch){
Tind <- c(1:V)[-chnind[,ch][!is.na(chnind[,ch])]]
if(length(Tind)==0){Tind <- c(1:V)}
chnind[v,ch] <- which(max(cor(truths[1:cctfit2$m, 1, v],truths[1:cctfit2$m, ch, Tind]))==cor(truths[1:cctfit2$m, 1, v],truths[1:cctfit2$m, ch, ]))
}}
}
nsamp <- cctfit$BUGSoutput$n.keep
if(cctfit$itemdiff==TRUE){
inds2 <- cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="lam")]]
inds2 <- matrix(inds2,cctfit2$m,cctfit2$V)
inds <- rbind(inds,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Tmu")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Ttau")]],
inds2,
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="lammu")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="lamtau")]],
matrix(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="gam")]],cctfit2$C-1,cctfit2$V),
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Emu")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Etau")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="amu")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="atau")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="bmu")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="btau")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="pi")]])
rm(inds2)
}else{
inds <- rbind(inds,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Tmu")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Ttau")]],
matrix(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="gam")]],cctfit2$C-1,cctfit2$V),
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Emu")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Etau")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="amu")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="atau")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="bmu")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="btau")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="pi")]])
}
einds <- cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]]
tmpeinds <- cctfit$BUGSoutput$sims.array[ ,,einds]
for(v in 1:cctfit$V){
for(ch in chstart:nch){
cctfit2$BUGSoutput$sims.array[ ,ch,inds[,v]] <- cctfit$BUGSoutput$sims.array[ ,ch,inds[,chnind[v,ch]]] #this cctfit2 vs. cctfit difference is intentional
if(chstart==2){cctfit2$BUGSoutput$sims.array[ ,ch,einds][cctfit$BUGSoutput$sims.array[ ,ch,einds]==chnind[v,ch]] <- chnind[v,1]}
if(chstart==1){cctfit2$BUGSoutput$sims.array[ ,ch,einds][cctfit$BUGSoutput$sims.array[ ,ch,einds]==v] <- chnind[v,1]}
}}
}
cctfit2$BUGSoutput$sims.list[["Om"]] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]]], c(nsamp*nch,cctfit$n))
cctfit2$BUGSoutput$mean$Om <- apply(cctfit2$BUGSoutput$sims.list[["Om"]],2,mean)
cctfit2$BUGSoutput$sims.matrix <- array(cctfit2$BUGSoutput$sims.array,c(nsamp*nch,dim(cctfit$BUGSoutput$sims.array)[3])) #this cctfit2 vs. cctfit difference is intentional
cctfit2$BUGSoutput$sims.list[["E"]] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="E")]]], c(nsamp*nch,cctfit$n))
cctfit2$BUGSoutput$sims.list[["a"]] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="a")]]], c(nsamp*nch,cctfit$n))
cctfit2$BUGSoutput$sims.list[["b"]] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="b")]]], c(nsamp*nch,cctfit$n))
#new code
cctfit2$BUGSoutput$mean[["E"]] <- apply(cctfit2$BUGSoutput$sims.list[["E"]],2,mean)
cctfit2$BUGSoutput$mean[["a"]] <- apply(cctfit2$BUGSoutput$sims.list[["a"]],2,mean)
cctfit2$BUGSoutput$mean[["b"]] <- apply(cctfit2$BUGSoutput$sims.list[["b"]],2,mean)
if(cctfit$itemdiff==TRUE){cctfit2$BUGSoutput$sims.list[["lam"]] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="lam")]]], c(nsamp*nch,cctfit$m))}
cctfit2$BUGSoutput$sims.list[["T"]] <- array(NA, c(nsamp*nch,cctfit$m,cctfit$V))
cctfit2$BUGSoutput$sims.list[["Tmu"]] <- array(NA, c(nsamp*nch,cctfit$V))
cctfit2$BUGSoutput$sims.list[["Ttau"]] <- array(NA, c(nsamp*nch,cctfit$V))
if(cctfit$itemdiff==TRUE){
cctfit2$BUGSoutput$sims.list[["lam"]] <- array(NA, c(nsamp*nch,cctfit$m,cctfit$V))
cctfit2$BUGSoutput$sims.list[["lammu"]] <- array(NA, c(nsamp*nch,cctfit$V))
cctfit2$BUGSoutput$sims.list[["lamtau"]] <- array(NA, c(nsamp*nch,cctfit$V))
}
cctfit2$BUGSoutput$sims.list[["gam"]] <- array(NA, c(nsamp*nch,cctfit$C-1,cctfit$V))
cctfit2$BUGSoutput$sims.list[["Emu"]] <- array(NA, c(nsamp*nch,cctfit$V))
cctfit2$BUGSoutput$sims.list[["Etau"]] <- array(NA, c(nsamp*nch,cctfit$V))
cctfit2$BUGSoutput$sims.list[["amu"]] <- array(NA, c(nsamp*nch,cctfit$V))
cctfit2$BUGSoutput$sims.list[["atau"]] <- array(NA, c(nsamp*nch,cctfit$V))
cctfit2$BUGSoutput$sims.list[["bmu"]] <- array(NA, c(nsamp*nch,cctfit$V))
cctfit2$BUGSoutput$sims.list[["btau"]] <- array(NA, c(nsamp*nch,cctfit$V))
cctfit2$BUGSoutput$sims.list[["pi"]] <- array(NA, c(nsamp*nch,cctfit$V))
if(cctfit$C==2 && length(dim(cctfit2$BUGSoutput$sims.list[["gam"]])) < 3 ){
cctfit2$BUGSoutput$sims.list[["gam"]] <- array(cctfit2$BUGSoutput$sims.list[["gam"]], c(nsamp*nch,cctfit$C-1,cctfit$V))
}
for(v in 1:cctfit2$V){
cctfit2$BUGSoutput$sims.list[["T"]][,,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="T")]][1:cctfit$m +((v-1)*cctfit$m)]], c(nsamp*nch,cctfit$m))
cctfit2$BUGSoutput$sims.list[["Tmu"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Tmu")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$sims.list[["Ttau"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Ttau")]][v]], c(nsamp*nch))
if(cctfit$itemdiff==TRUE){
cctfit2$BUGSoutput$sims.list[["lam"]][,,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="lam")]][1:cctfit$m +((v-1)*cctfit$m)]], c(nsamp*nch,cctfit$m))
cctfit2$BUGSoutput$sims.list[["lammu"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="lammu")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$sims.list[["lamtau"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="lamtau")]][v]], c(nsamp*nch))
}
cctfit2$BUGSoutput$sims.list[["gam"]][,,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="gam")]][1:(cctfit$C-1) +((v-1)*(cctfit$C-1))]], c(nsamp*nch,cctfit$C-1))
cctfit2$BUGSoutput$sims.list[["Emu"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Emu")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$sims.list[["Etau"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Etau")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$sims.list[["amu"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="amu")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$sims.list[["atau"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="atau")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$sims.list[["bmu"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="bmu")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$sims.list[["btau"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="btau")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$sims.list[["pi"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="pi")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$mean$T[,v] <- apply(cctfit2$BUGSoutput$sims.list[["T"]][,,v],2,mean)
cctfit2$BUGSoutput$mean$Tmu[v] <- mean(cctfit2$BUGSoutput$sims.list[["Tmu"]][,v])
cctfit2$BUGSoutput$mean$Ttau[v] <- mean(cctfit2$BUGSoutput$sims.list[["Ttau"]][,v])
if(cctfit$itemdiff==TRUE){
cctfit2$BUGSoutput$mean$lam[,v] <- apply(cctfit2$BUGSoutput$sims.list[["lam"]][,,v],2,mean)
cctfit2$BUGSoutput$mean$lammu[v] <- mean(cctfit2$BUGSoutput$sims.list[["lammu"]][,v])
cctfit2$BUGSoutput$mean$lamtau[v] <- mean(cctfit2$BUGSoutput$sims.list[["lamtau"]][,v])
}
if(cctfit$C==2){
cctfit2$BUGSoutput$mean$gam[v] <- mean(cctfit2$BUGSoutput$sims.list[["gam"]][,,v])
}else{cctfit2$BUGSoutput$mean$gam[,v] <- apply(cctfit2$BUGSoutput$sims.list[["gam"]][,,v],2,mean)}
cctfit2$BUGSoutput$mean$Emu[v] <- mean(cctfit2$BUGSoutput$sims.list[["Emu"]][,v])
cctfit2$BUGSoutput$mean$Etau[v] <- mean(cctfit2$BUGSoutput$sims.list[["Etau"]][,v])
cctfit2$BUGSoutput$mean$amu[v] <- mean(cctfit2$BUGSoutput$sims.list[["amu"]][,v])
cctfit2$BUGSoutput$mean$atau[v] <- mean(cctfit2$BUGSoutput$sims.list[["atau"]][,v])
cctfit2$BUGSoutput$mean$bmu[v] <- mean(cctfit2$BUGSoutput$sims.list[["bmu"]][,v])
cctfit2$BUGSoutput$mean$btau[v] <- mean(cctfit2$BUGSoutput$sims.list[["btau"]][,v])
cctfit2$BUGSoutput$mean$pi[v] <- mean(cctfit2$BUGSoutput$sims.list[["pi"]][,v])
}
if(nch!=1){
cctfit2$BUGSoutput$summary[,1] <- apply(apply(cctfit2$BUGSoutput$sims.array,c(2,3),mean),2,mean)
cctfit2$BUGSoutput$summary[,2] <- apply(apply(cctfit2$BUGSoutput$sims.array,c(2,3),sd),2,sd)
cctfit2$BUGSoutput$summary[,3:7] <- t(apply(cctfit2$BUGSoutput$sims.matrix,2,function(x) quantile(x,probs=c(.025,.25,.50,.75,.975))))
cctfit2$BUGSoutput$summary[,8] <- Rhat(cctfit2$BUGSoutput$sims.array)
cctfit2$BUGSoutput$summary[,8][is.nan(cctfit2$BUGSoutput$summary[,8])] <- 1.000000
dimnames(cctfit2$BUGSoutput$sims.array) <- list(NULL,NULL,rownames(cctfit$BUGSoutput$summary))
dimnames(cctfit2$BUGSoutput$sims.matrix) <- list(NULL,rownames(cctfit$BUGSoutput$summary))
}else{
cctfit2$BUGSoutput$summary[,1] <- apply(cctfit2$BUGSoutput$sims.array,2,mean)
cctfit2$BUGSoutput$summary[,2] <- apply(cctfit2$BUGSoutput$sims.array,2,sd)
cctfit2$BUGSoutput$summary[,3:7] <- t(apply(cctfit2$BUGSoutput$sims.matrix,2,function(x) quantile(x,probs=c(.025,.25,.50,.75,.975))))
cctfit2$BUGSoutput$summary <- cctfit2$BUGSoutput$summary[,-c(8,9)]
dimnames(cctfit2$BUGSoutput$sims.array) <- list(NULL,NULL,rownames(cctfit$BUGSoutput$summary))
dimnames(cctfit2$BUGSoutput$sims.matrix) <- list(NULL,rownames(cctfit$BUGSoutput$summary))
}
cctfit <- cctfit2; rm(cctfit2)
return(cctfit)
}
######################
#Algorithm that corrects for label switching for the CRM
######################
labelswitchalgcrm <- function(cctfit,chnind=0){
cctfit2 <- cctfit
nch <- cctfit2$BUGSoutput$n.chains
nsamp <- cctfit2$BUGSoutput$n.keep
if(nch!=1){
ntruths <- cctfit2$V
truths <- array(NA,c(cctfit2$m,nch,ntruths))
inds <- cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="T")]]
inds <- matrix(inds,cctfit2$m,cctfit2$V)
for(v in 1:cctfit$V){truths[,,v] <- t(apply(cctfit$BUGSoutput$sims.array[ ,,inds[,v]],c(2,3),mean))}
V <- cctfit2$V
chstart <- 1
if(length(chnind)==1){
chstart <- 2
chnind <- array(NA,c(V,nch))
chnind[1:cctfit2$V,1] <- 1:cctfit2$V
for(v in 1:cctfit2$V){
for(ch in chstart:nch){
Tind <- c(1:cctfit2$V)[-chnind[,ch][!is.na(chnind[,ch])]]
if(length(Tind)==0){Tind <- c(1:cctfit2$V)}
chnind[v,ch] <- which(max(cor(truths[1:cctfit2$m, 1, v],truths[1:cctfit2$m, ch, Tind]))==cor(truths[1:cctfit2$m, 1, v],truths[1:cctfit2$m, ch, ]))
}}
}
nsamp <- cctfit$BUGSoutput$n.keep
if(cctfit$itemdiff==TRUE){
inds2 <- cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="lam")]]
inds2 <- matrix(inds2,cctfit2$m,cctfit2$V)
inds <- rbind(inds,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Tmu")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Ttau")]],
inds2,
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="lammu")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="lamtau")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Emu")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Etau")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="amu")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="atau")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="bmu")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="btau")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="pi")]])
rm(inds2)
}else{
inds <- rbind(inds,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Tmu")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Ttau")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Emu")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Etau")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="amu")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="atau")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="bmu")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="btau")]],
cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="pi")]])
}
einds <- cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]]
tmpeinds <- cctfit$BUGSoutput$sims.array[ ,,einds]
for(v in 1:cctfit$V){
for(ch in chstart:nch){
cctfit2$BUGSoutput$sims.array[ ,ch,inds[,v]] <- cctfit$BUGSoutput$sims.array[ ,ch,inds[,chnind[v,ch]]] #this cctfit2 vs. cctfit difference is intentional
if(chstart==2){cctfit2$BUGSoutput$sims.array[ ,ch,einds][cctfit$BUGSoutput$sims.array[ ,ch,einds]==chnind[v,ch]] <- chnind[v,1]}
if(chstart==1){cctfit2$BUGSoutput$sims.array[ ,ch,einds][cctfit$BUGSoutput$sims.array[ ,ch,einds]==v] <- chnind[v,1]}
}}
}
cctfit2$BUGSoutput$sims.list[["Om"]] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]]], c(nsamp*nch,cctfit$n))
cctfit2$BUGSoutput$mean$Om <- apply(cctfit2$BUGSoutput$sims.list[["Om"]],2,mean)
cctfit2$BUGSoutput$sims.matrix <- array(cctfit2$BUGSoutput$sims.array,c(nsamp*nch,dim(cctfit$BUGSoutput$sims.array)[3])) #this cctfit2 vs. cctfit difference is intentional
cctfit2$BUGSoutput$sims.list[["E"]] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="E")]]], c(nsamp*nch,cctfit$n))
cctfit2$BUGSoutput$sims.list[["a"]] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="a")]]], c(nsamp*nch,cctfit$n))
cctfit2$BUGSoutput$sims.list[["b"]] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="b")]]], c(nsamp*nch,cctfit$n))
#new code
cctfit2$BUGSoutput$mean[["E"]] <- apply(cctfit2$BUGSoutput$sims.list[["E"]],2,mean)
cctfit2$BUGSoutput$mean[["a"]] <- apply(cctfit2$BUGSoutput$sims.list[["a"]],2,mean)
cctfit2$BUGSoutput$mean[["b"]] <- apply(cctfit2$BUGSoutput$sims.list[["b"]],2,mean)
if(cctfit$itemdiff==TRUE){cctfit2$BUGSoutput$sims.list[["lam"]] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="lam")]]], c(nsamp*nch,cctfit$m))}
cctfit2$BUGSoutput$sims.list[["T"]] <- array(NA, c(nsamp*nch,cctfit$m,cctfit$V))
cctfit2$BUGSoutput$sims.list[["Tmu"]] <- array(NA, c(nsamp*nch,cctfit$V))
cctfit2$BUGSoutput$sims.list[["Ttau"]] <- array(NA, c(nsamp*nch,cctfit$V))
if(cctfit$itemdiff==TRUE){
cctfit2$BUGSoutput$sims.list[["lam"]] <- array(NA, c(nsamp*nch,cctfit$m,cctfit$V))
cctfit2$BUGSoutput$sims.list[["lammu"]] <- array(NA, c(nsamp*nch,cctfit$V))
cctfit2$BUGSoutput$sims.list[["lamtau"]] <- array(NA, c(nsamp*nch,cctfit$V))
}
cctfit2$BUGSoutput$sims.list[["Emu"]] <- array(NA, c(nsamp*nch,cctfit$V))
cctfit2$BUGSoutput$sims.list[["Etau"]] <- array(NA, c(nsamp*nch,cctfit$V))
cctfit2$BUGSoutput$sims.list[["amu"]] <- array(NA, c(nsamp*nch,cctfit$V))
cctfit2$BUGSoutput$sims.list[["atau"]] <- array(NA, c(nsamp*nch,cctfit$V))
cctfit2$BUGSoutput$sims.list[["bmu"]] <- array(NA, c(nsamp*nch,cctfit$V))
cctfit2$BUGSoutput$sims.list[["btau"]] <- array(NA, c(nsamp*nch,cctfit$V))
cctfit2$BUGSoutput$sims.list[["pi"]] <- array(NA, c(nsamp*nch,cctfit$V))
for(v in 1:cctfit2$V){
cctfit2$BUGSoutput$sims.list[["T"]][,,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="T")]][1:cctfit$m +((v-1)*cctfit$m)]], c(nsamp*nch,cctfit$m))
cctfit2$BUGSoutput$sims.list[["Tmu"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Tmu")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$sims.list[["Ttau"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Ttau")]][v]], c(nsamp*nch))
if(cctfit$itemdiff==TRUE){
cctfit2$BUGSoutput$sims.list[["lam"]][,,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="lam")]][1:cctfit$m +((v-1)*cctfit$m)]], c(nsamp*nch,cctfit$m))
cctfit2$BUGSoutput$sims.list[["lammu"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="lammu")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$sims.list[["lamtau"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="lamtau")]][v]], c(nsamp*nch))
}
cctfit2$BUGSoutput$sims.list[["Emu"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Emu")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$sims.list[["Etau"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Etau")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$sims.list[["amu"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="amu")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$sims.list[["atau"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="atau")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$sims.list[["bmu"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="bmu")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$sims.list[["btau"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="btau")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$sims.list[["pi"]][,v] <- array(cctfit2$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="pi")]][v]], c(nsamp*nch))
cctfit2$BUGSoutput$mean$T[,v] <- apply(cctfit2$BUGSoutput$sims.list[["T"]][,,v],2,mean)
cctfit2$BUGSoutput$mean$Tmu[v] <- mean(cctfit2$BUGSoutput$sims.list[["Tmu"]][,v])
cctfit2$BUGSoutput$mean$Ttau[v] <- mean(cctfit2$BUGSoutput$sims.list[["Ttau"]][,v])
if(cctfit$itemdiff==TRUE){
cctfit2$BUGSoutput$mean$lam[,v] <- apply(cctfit2$BUGSoutput$sims.list[["lam"]][,,v],2,mean)
cctfit2$BUGSoutput$mean$lammu[v] <- mean(cctfit2$BUGSoutput$sims.list[["lammu"]][,v])
cctfit2$BUGSoutput$mean$lamtau[v] <- mean(cctfit2$BUGSoutput$sims.list[["lamtau"]][,v])
}
cctfit2$BUGSoutput$mean$Emu[v] <- mean(cctfit2$BUGSoutput$sims.list[["Emu"]][,v])
cctfit2$BUGSoutput$mean$Etau[v] <- mean(cctfit2$BUGSoutput$sims.list[["Etau"]][,v])
cctfit2$BUGSoutput$mean$amu[v] <- mean(cctfit2$BUGSoutput$sims.list[["amu"]][,v])
cctfit2$BUGSoutput$mean$atau[v] <- mean(cctfit2$BUGSoutput$sims.list[["atau"]][,v])
cctfit2$BUGSoutput$mean$bmu[v] <- mean(cctfit2$BUGSoutput$sims.list[["bmu"]][,v])
cctfit2$BUGSoutput$mean$btau[v] <- mean(cctfit2$BUGSoutput$sims.list[["btau"]][,v])
cctfit2$BUGSoutput$mean$pi[v] <- mean(cctfit2$BUGSoutput$sims.list[["pi"]][,v])
}
if(nch!=1){
cctfit2$BUGSoutput$summary[,1] <- apply(apply(cctfit2$BUGSoutput$sims.array,c(2,3),mean),2,mean)
cctfit2$BUGSoutput$summary[,2] <- apply(apply(cctfit2$BUGSoutput$sims.array,c(2,3),sd),2,sd)
cctfit2$BUGSoutput$summary[,3:7] <- t(apply(cctfit2$BUGSoutput$sims.matrix,2,function(x) quantile(x,probs=c(.025,.25,.50,.75,.975))))
cctfit2$BUGSoutput$summary[,8] <- Rhat(cctfit2$BUGSoutput$sims.array)
cctfit2$BUGSoutput$summary[,8][is.nan(cctfit2$BUGSoutput$summary[,8])] <- 1.000000
dimnames(cctfit2$BUGSoutput$sims.array) <- list(NULL,NULL,rownames(cctfit$BUGSoutput$summary))
dimnames(cctfit2$BUGSoutput$sims.matrix) <- list(NULL,rownames(cctfit$BUGSoutput$summary))
}else{
cctfit2$BUGSoutput$summary[,1] <- apply(cctfit2$BUGSoutput$sims.array,2,mean)
cctfit2$BUGSoutput$summary[,2] <- apply(cctfit2$BUGSoutput$sims.array,2,sd)
cctfit2$BUGSoutput$summary[,3:7] <- t(apply(cctfit2$BUGSoutput$sims.matrix,2,function(x) quantile(x,probs=c(.025,.25,.50,.75,.975))))
cctfit2$BUGSoutput$summary <- cctfit2$BUGSoutput$summary[,-c(8,9)]
dimnames(cctfit2$BUGSoutput$sims.array) <- list(NULL,NULL,rownames(cctfit$BUGSoutput$summary))
dimnames(cctfit2$BUGSoutput$sims.matrix) <- list(NULL,rownames(cctfit$BUGSoutput$summary))
}
cctfit <- cctfit2; rm(cctfit2)
return(cctfit)
}
######################
#This is run after the jags inference, we are still in the 'Apply CCT Model' function
#- Reports the number of Rhats above 1.05 and 1.10 (before applying the algorithm that corrects for label-switching)
#- Detects if a mixture model was applied
#- If so, the appropriate label-correcting algorithm is applied
#- Recalculates the Rhats, DIC, and statistics for the parameters
#- Reports the number of Rhats above 1.05 and 1.10
#- Outputs the DIC that is calculated after the label-correcting algorithm (if applicable)
######################
# message("\n ...Inference complete, data is saved as '",guidat$varname,"'")
# if(clusters > 1){
# message("\n 'cctfit$respmem' provides the respondent clustering")
# }
hdi <- function(sampleVec, credMass=0.95){
sortedPts = sort(sampleVec)
ciIdxInc = floor(credMass * length( sortedPts) )
nCIs = length( sortedPts) - ciIdxInc
ciwidth = rep(0, nCIs)
for(i in 1:nCIs){
ciwidth[i] = sortedPts[i + ciIdxInc] - sortedPts[i]}
HDImin = sortedPts[which.min(ciwidth)]
HDImax = sortedPts[which.min(ciwidth)+ciIdxInc]
HDIlim = c(HDImin,HDImax)
return(HDIlim)
}
message("\n ...Performing final calculations")
# }
if(cctfit$BUGSoutput$n.chains > 1){
cctfit$BUGSoutput$summary[,8] <- Rhat(cctfit$BUGSoutput$sims.array)
cctfit$BUGSoutput$summary[,8][is.nan(cctfit$BUGSoutput$summary[,8])] <- 1.000000
if(cctfit$whmodel== "GCM"){
cctfit$Rhat$ncp <- length(cctfit$BUGSoutput$summary[-c(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]],cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Z")]]),8])
cctfit$Rhat$above110 <- sum(cctfit$BUGSoutput$summary[-c(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]],cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Z")]]),8]>1.10)
cctfit$Rhat$above105 <- sum(cctfit$BUGSoutput$summary[-c(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]],cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Z")]]),8]>1.050)
}else{
cctfit$Rhat$ncp <- length(cctfit$BUGSoutput$summary[-c(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]]),8])
cctfit$Rhat$above110 <- sum(cctfit$BUGSoutput$summary[-c(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]]),8]>1.10)
cctfit$Rhat$above105 <- sum(cctfit$BUGSoutput$summary[-c(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]]),8]>1.050)
}
message(paste("\nFor Continuous Parameters"))
message(paste("Number of Rhats above 1.10 : ",cctfit$Rhat$above110,"/",cctfit$Rhat$ncp,"\nNumber of Rhats above 1.05 : ",cctfit$Rhat$above105,"/",cctfit$Rhat$ncp))
}
if(cctfit$V>1 && cctfit$BUGSoutput$n.chains==1){
tmp <- unique(apply(cctfit$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]]],c(2),Mode))
if(length(tmp)==1){
message(paste("\n ...This chain has ",length(tmp)," culture rather than the ",cctfit$V," cultures requested", sep=""))
message(paste("\n ...Try running the inference again",sep="" ))
}
if(length(tmp)!=1 && length(tmp) < cctfit$V){
message(paste("\n ...This chain has ",tmp," cultures rather than the ",cctfit$V," cultures requested", sep=""))
message(paste("\n ...Try running the inference again",sep="" ))
}
}
if(cctfit$V>1 && cctfit$BUGSoutput$n.chains > 1){
message("\n ...More than 1 culture applied with more than 1 chain")
einds <- cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]]
tmp <- apply(apply(cctfit$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]]],c(2,3),Mode),1,unique)
chntoremove <- NULL
if(is.list(tmp)){
for(i in 1:cctfit$BUGSoutput$n.chains){
if(length(tmp[[i]])!=cctfit$V){chntoremove <- c(chntoremove,i)}
}
}
if(length(dim(tmp))==2){
for(i in 1:cctfit$BUGSoutput$n.chains){
if(length(tmp[,i])!=cctfit$V){chntoremove <- c(chntoremove,i)}
}
}
if(is.null(dim(tmp)) && !is.list(tmp)){chntoremove <- c(1:cctfit$BUGSoutput$n.chains)}
if(length(chntoremove)>0){
if(length(chntoremove) < cctfit$BUGSoutput$n.chains){
if(length(chntoremove)==1){
message(paste("\n ...",length(chntoremove)," chain out of ",cctfit$BUGSoutput$n.chains," had fewer than ",cctfit$V," cultures requested",sep=""))
message(paste("\n ...", "removing the ", length(chntoremove)," chain", sep=""))
}else{
message(paste("\n ...",length(chntoremove)," chains out of ",cctfit$BUGSoutput$n.chains," had fewer than ",cctfit$V," cultures requested",sep=""))
message(paste("\n ...", "removing these ", length(chntoremove)," chains", sep=""))
}
cctfit$BUGSoutput$n.chains <- cctfit$BUGSoutput$n.chains-length(chntoremove)
cctfit$BUGSoutput$n.chain <- cctfit$BUGSoutput$n.chains
cctfit$BUGSoutput$n.sims <- cctfit$BUGSoutput$n.chains*cctfit$BUGSoutput$n.keep
if(cctfit$BUGSoutput$n.chain==1){
cctfit$BUGSoutput$sims.array <- array(cctfit$BUGSoutput$sims.array[,-chntoremove,], c(dim(cctfit$BUGSoutput$sims.array)[1],1,dim(cctfit$BUGSoutput$sims.array)[3]))
}else{cctfit$BUGSoutput$sims.array <- cctfit$BUGSoutput$sims.array[,-chntoremove,]}
}else{
message(paste("\n ...All chains out of ",cctfit$BUGSoutput$n.chains," had fewer than ",cctfit$V," cultures requested", sep=""))
message(paste("\n ...Try running the inference again",sep="" ))
}
}
message("\n ...Computing the most-consistent labeling across chains")
if(cctfit$whmodel=="GCM"){
cctfit <- labelswitchalggcm(cctfit)
cctfit$respmem <- apply(cctfit$BUGSoutput$sims.list$Om[,],2,Mode)
tmeans <- cctfit$BUGSoutput$mean$Z
tmeans[tmeans<.5] <- tmeans[tmeans<.5]+1
ind <- rank(apply(abs(1-tmeans),2,mean))
}
if(cctfit$whmodel=="LTRM"){
cctfit <- labelswitchalgltrm(cctfit)
cctfit$respmem <- apply(cctfit$BUGSoutput$sims.list$Om[,],2,Mode)
tmeans <- cctfit$BUGSoutput$mean$T
ind <- rank(-apply(tmeans,2,sd))
}
if(cctfit$whmodel=="CRM"){
cctfit <- labelswitchalgcrm(cctfit)
cctfit$respmem <- apply(cctfit$BUGSoutput$sims.list$Om[,],2,Mode)
tmeans <- cctfit$BUGSoutput$mean$T
ind <- rank(-apply(tmeans,2,sd))
}
# old code:
# if(cctfit$BUGSoutput$n.chains > 1){
# message(paste("\nFor Continuous Parameters:"))
# if(cctfit$whmodel== "GCM"){
# message(paste("Number of Rhats above 1.10 : ",sum(cctfit$BUGSoutput$summary[-c(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]],cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Z")]]),8]>1.10),"/",length(cctfit$BUGSoutput$summary[-c(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]],cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Z")]]),8]),"\nNumber of Rhats above 1.05 : ",sum(cctfit$BUGSoutput$summary[-c(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]],cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Z")]]),8]>1.050),"/",length(cctfit$BUGSoutput$summary[-c(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]],cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Z")]]),8]),sep=""))
# }else{
# message(paste("Number of Rhats above 1.10 : ",sum(cctfit$BUGSoutput$summary[-c(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]]),8]>1.10),"/",length(cctfit$BUGSoutput$summary[-c(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]]),8]),"\nNumber of Rhats above 1.05 : ",sum(cctfit$BUGSoutput$summary[-c(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]]),8]>1.050),"/",length(cctfit$BUGSoutput$summary[-c(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]]),8]),sep=""))
# }
# }
if(cctfit$BUGSoutput$n.chains > 1){
cctfit$BUGSoutput$summary[,8] <- Rhat(cctfit$BUGSoutput$sims.array)
cctfit$BUGSoutput$summary[,8][is.nan(cctfit$BUGSoutput$summary[,8])] <- 1.000000
if(cctfit$whmodel== "GCM"){
cctfit$Rhat$ncp <- length(cctfit$BUGSoutput$summary[-c(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]],cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Z")]]),8])
cctfit$Rhat$above110 <- sum(cctfit$BUGSoutput$summary[-c(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]],cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Z")]]),8]>1.10)
cctfit$Rhat$above105 <- sum(cctfit$BUGSoutput$summary[-c(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]],cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Z")]]),8]>1.050)
}else{
cctfit$Rhat$ncp <- length(cctfit$BUGSoutput$summary[-c(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]]),8])
cctfit$Rhat$above110 <- sum(cctfit$BUGSoutput$summary[-c(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]]),8]>1.10)
cctfit$Rhat$above105 <- sum(cctfit$BUGSoutput$summary[-c(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]]),8]>1.050)
}
message(paste("\nFor Continuous Parameters"))
message(paste("Number of Rhats above 1.10 : ",cctfit$Rhat$above110,"/",cctfit$Rhat$ncp,"\nNumber of Rhats above 1.05 : ",cctfit$Rhat$above105,"/",cctfit$Rhat$ncp))
}
if(cctfit$BUGSoutput$n.chains > 1){
if(gui==TRUE){message(paste("\nFor Discrete Parameters:"))
message(paste("Use button 'Traceplot Discrete' to see their trace plots",sep=""))
}else{
message(paste("\nFor Discrete Parameters:"))
message(paste("Use function 'dtraceplot()' to see their trace plots",sep=""))
message(paste("\nUse function 'cctmemb()' to see the respondent cluster assignments",sep=""))
}
}
}else{
cctfit$respmem <- rep(1,cctfit$n)
if(cctfit$BUGSoutput$n.chains > 1){
if(gui==TRUE){message(paste("\nFor Discrete Parameters:"))
message(paste("Use button 'Traceplot Discrete' to see their trace plots",sep=""))
}else{
message(paste("\nFor Discrete Parameters:"))
message(paste("Use function 'dtraceplot()' to see their trace plots",sep=""))
}
}
}
######################
#Calculates Read-out tables for subjects and items
#
#######################
Mode <- function(x) {ux <- unique(x); ux[which.max(tabulate(match(x, ux)))] }
cctfit$paramsall <- cctfit$parameters.to.save
cctfit$params <- cctfit$paramsall[cctfit$parameters.to.save %in% c("Z","T","Om","th","E","a","b","g","lam","gam")]
if(cctfit$V > 1){appdims <- c(2,3)}else{appdims <- 2}
a <- data.frame(ans=round(cctfit$BUGSoutput$mean[[cctfit$params[cctfit$params %in% c("Z","T")]]],2));
#a <- data.frame(ans=round(apply(cctfit$BUGSoutput$sims.list[[cctfit$params[1]]],c(2,3),mean),2)) # should be c(2,3) in all cases here;
dimnames(a)[[2]] <- gsub(pattern="s.",replacement=paste("s",sep=""),x=dimnames(a)[[2]]); dimnames(a)[[2]] <- paste(dimnames(a)[[2]],"_",cctfit$params[cctfit$params %in% c("Z","T")],sep="")
a <- data.frame(item=1:cctfit$m,a)
if(cctfit$itemdiff==TRUE){
d <- data.frame(diff=round(cctfit$BUGSoutput$mean[["lam"]],2))
#d <- data.frame(diff=round(apply(cctfit$BUGSoutput$sims.list[["lam"]],c(2,3),mean),2)); # should be c(2,3) in all cases here;
dimnames(d)[[2]] <- gsub(pattern="ff.",replacement=paste("ff",sep=""),x=dimnames(d)[[2]]); dimnames(d)[[2]] <- paste(dimnames(d)[[2]],"_lam",sep="")
a <- data.frame(a,d)
}
cctfit$item <- a
p <- data.frame(group_Om=apply(cctfit$BUGSoutput$sims.list$Om[,],2,Mode),
#comp=apply(cctfit$BUGSoutput$sims.list[[cctfit$params[3]]],2,mean));
comp=round(cctfit$BUGSoutput$mean[[cctfit$params[cctfit$params %in% c("th","E")]]],2));
dimnames(p)[[2]][2] <- paste(dimnames(p)[[2]][2],"_",cctfit$params[cctfit$params %in% c("th","E")],sep="")
b <- round(as.data.frame(cctfit$BUGSoutput$mean[which(names(cctfit$BUGSoutput$mean) %in% c("a","b","g"))]),2)
dimnames(b)[[2]] <- paste("bias_",dimnames(b)[[2]],sep="")
p <- data.frame(participant=1:cctfit$n,p,b)
cctfit$subj <- p
tmphdi <- apply(cctfit$BUGSoutput$sims.list[[cctfit$params[1]]],appdims,hdi)
ahdi <- NULL;
if(cctfit$V>1){for(i in 1:cctfit$V){ ahdi <- cbind(ahdi,t(tmphdi[,,i])) };
}else{ahdi <- t(tmphdi)}; ahdi <- data.frame(ahdi)
for(i in seq(from=1,to=(cctfit$V*2),by=2)){dimnames(ahdi)[[2]][i] <- paste("l_ans",ceiling(i/2),sep=""); dimnames(ahdi)[[2]][i+1] <- paste("u_ans",ceiling(i/2),sep="")}
if(cctfit$itemdiff==TRUE){
tmphdi <- apply(cctfit$BUGSoutput$sims.list[["lam"]],appdims,hdi)
dhdi <- NULL;
if(cctfit$V>1){
for(i in 1:cctfit$V){ dhdi <- round(cbind(dhdi,t(tmphdi[,,i])),2) };
}else{dhdi <- t(tmphdi)}; dhdi <- data.frame(dhdi)
for(i in seq(from=1,to=(cctfit$V*2),by=2)){dimnames(dhdi)[[2]][i] <- paste("l_diff",ceiling(i/2),sep=""); dimnames(dhdi)[[2]][i+1] <- paste("u_diff",ceiling(i/2),sep="")}
ahdi <- cbind(ahdi,dhdi)
}
cctfit$itemhdi <- data.frame(item=1:cctfit$m,ahdi);
tmphdi <- NULL; tmpn <- NULL; pp <- cctfit$params[cctfit$params %in% c("th","E","a","b","g")]
for(i in 1:length(pp)){
tmphdi <- cbind(tmphdi,round(t(apply(cctfit$BUGSoutput$sims.list[[pp[i]]],2,hdi)),2))
tmpn <- c(tmpn,c(paste("l_",pp[i],sep=""),paste("u_",pp[i],sep="")))
}
tmphdi <- data.frame(tmphdi); dimnames(tmphdi)[[2]] <- tmpn;
cctfit$subjhdi <- data.frame(participant=1:cctfit$n,tmphdi)
if(cctfit$whmodel=="LTRM"){
pp <- "gam"
tmphdi <- apply(cctfit$BUGSoutput$sims.list[[pp]],appdims,hdi)
ahdi <- NULL;
a <- data.frame(gam=round(cctfit$BUGSoutput$mean[[pp]],2));
dimnames(a)[[2]] <- gsub(pattern="m.",replacement=paste("m",sep=""),x=dimnames(a)[[2]]);
a <- data.frame(boundary=1:(cctfit$C-1),a)
tmphdi <- apply(cctfit$BUGSoutput$sims.list[[pp]],appdims,hdi)
ahdi <- NULL;
if(cctfit$V>1){for(i in 1:cctfit$V){ ahdi <- cbind(ahdi,t(tmphdi[,,i])) };
}else{ahdi <- t(tmphdi)}; ahdi <- data.frame(ahdi)
for(i in seq(from=1,to=(cctfit$V*2),by=2)){dimnames(ahdi)[[2]][i] <- paste("l_gam",ceiling(i/2),sep=""); dimnames(ahdi)[[2]][i+1] <- paste("u_gam",ceiling(i/2),sep="")}
a <- data.frame(a,round(ahdi,digits=2))
cctfit$cat <- a
}
######################
#DIC is calculated for the model that was applied
#The likelihood of the model is evaluated at each node of each sample and saved to cctfit$Lik
######################
message("\n ...Calculating DIC")
cctfit$BUGSoutput$pD <- var(cctfit$BUGSoutput$sims.list$deviance[,1])/2
nsimstouse <- min(cctfit$BUGSoutput$n.sims,1000)
ind <- sample(1:cctfit$BUGSoutput$n.sims,nsimstouse)
storelik <- 0
if(cctfit$whmodel=="GCM"){
if(cctfit$V==1){
cctfit$V <- 1;
cctfit$BUGSoutput$sims.list$Om <- array(1,c(cctfit$BUGSoutput$n.sims,cctfit$n));
if(length(dim(cctfit$BUGSoutput$sims.list$Z)) < 3){
cctfit$BUGSoutput$sims.list$Z <- array(cctfit$BUGSoutput$sims.list[["Z"]][,],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))}
if(cctfit$itemdiff==TRUE){
if(length(dim(cctfit$BUGSoutput$sims.list$lam)) < 3){
cctfit$BUGSoutput$sims.list$lam <- array(cctfit$BUGSoutput$sims.list[["lam"]][,],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))}
}
}
if(cctfit$itemdiff==FALSE){
if(cctfit$V==1){
cctfit$BUGSoutput$sims.list$lam <- array(.5,c(cctfit$BUGSoutput$n.sims,cctfit$m,1))
}
if(cctfit$V > 1){
cctfit$BUGSoutput$sims.list$lam <- array(.5,c(cctfit$BUGSoutput$n.sims,cctfit$m,cctfit$V))
}
}
cctfit$BUGSoutput$sims.list$D <- array(-1,c(cctfit$BUGSoutput$n.sims,cctfit$n,cctfit$m))
storefulllik <- 0
if(storefulllik==1){
cctfit$Lik <- array(NA, c(cctfit$n,cctfit$m,cctfit$BUGSoutput$n.sims));
for(samp in 1:cctfit$BUGSoutput$n.sims){
cctfit$BUGSoutput$sims.list[["D"]][samp,,] <- (cctfit$BUGSoutput$sims.list[["th"]][samp,]*t(1-cctfit$BUGSoutput$sims.list[["lam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])) /
( (cctfit$BUGSoutput$sims.list[["th"]][samp,]*t(1-cctfit$BUGSoutput$sims.list[["lam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])) +
((1-cctfit$BUGSoutput$sims.list[["th"]][samp,])*t(cctfit$BUGSoutput$sims.list[["lam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])) )
cctfit$Lik[,,samp] <- t( ( (t(cctfit$BUGSoutput$sims.list[["D"]][samp,,])+(t(1-cctfit$BUGSoutput$sims.list[["D"]][samp,,])%*%diag(cctfit$BUGSoutput$sims.list[["g"]][samp,])))^(t(cctfit$data[,])*cctfit$BUGSoutput$sims.list[["Z"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]]))*(t(1-cctfit$BUGSoutput$sims.list[["D"]][samp,,])%*%diag(cctfit$BUGSoutput$sims.list[["g"]][samp,]))^(t(cctfit$data[,])*(1-cctfit$BUGSoutput$sims.list[["Z"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]]))*(t(1-cctfit$BUGSoutput$sims.list[["D"]][samp,,])%*%diag(1-cctfit$BUGSoutput$sims.list[["g"]][samp,]))^((1-t(cctfit$data[,]))*cctfit$BUGSoutput$sims.list[["Z"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])*(t(cctfit$BUGSoutput$sims.list[["D"]][samp,,])+t(1-cctfit$BUGSoutput$sims.list[["D"]][samp,,])%*%diag(1-cctfit$BUGSoutput$sims.list[["g"]][samp,]))^((1-t(cctfit$data[,]))*(1-cctfit$BUGSoutput$sims.list[["Z"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])) )
if(cctfit$mval==TRUE){cctfit$Lik[cbind(datob$thena,samp)] <- 1}
}
cctfit$BUGSoutput$sims.list$deviance <- array(-2*apply(log(cctfit$Lik),3,sum),c(cctfit$BUGSoutput$n.sims,1))
}else{
cctfit$LogLik <- array(NA, c(cctfit$BUGSoutput$n.sims));
for(samp in 1:cctfit$BUGSoutput$n.sims){
Dtmp <- (cctfit$BUGSoutput$sims.list[["th"]][samp,]*t(1-cctfit$BUGSoutput$sims.list[["lam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])) /
( (cctfit$BUGSoutput$sims.list[["th"]][samp,]*t(1-cctfit$BUGSoutput$sims.list[["lam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])) +
((1-cctfit$BUGSoutput$sims.list[["th"]][samp,])*t(cctfit$BUGSoutput$sims.list[["lam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])) )
Liktmp <- t( ( (t(Dtmp)+(t(1-Dtmp)%*%diag(cctfit$BUGSoutput$sims.list[["g"]][samp,])))^(t(cctfit$data[,])*cctfit$BUGSoutput$sims.list[["Z"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]]))*(t(1-Dtmp)%*%diag(cctfit$BUGSoutput$sims.list[["g"]][samp,]))^(t(cctfit$data[,])*(1-cctfit$BUGSoutput$sims.list[["Z"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]]))*(t(1-Dtmp)%*%diag(1-cctfit$BUGSoutput$sims.list[["g"]][samp,]))^((1-t(cctfit$data[,]))*cctfit$BUGSoutput$sims.list[["Z"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])*(t(Dtmp)+t(1-Dtmp)%*%diag(1-cctfit$BUGSoutput$sims.list[["g"]][samp,]))^((1-t(cctfit$data[,]))*(1-cctfit$BUGSoutput$sims.list[["Z"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])) )
if(cctfit$mval==TRUE){Liktmp[datob$thena] <- 1}
cctfit$LogLik[samp] <- sum(log(Liktmp))
}
cctfit$BUGSoutput$sims.list$deviance <- array(-2*(cctfit$LogLik),c(cctfit$BUGSoutput$n.sims,1))
}
cctfit$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="deviance")]]] <- array(cctfit$BUGSoutput$sims.list$deviance,c(cctfit$BUGSoutput$n.keep,cctfit$BUGSoutput$n.chains))
cctfit$BUGSoutput$sims.matrix[,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="deviance")]]] <- cctfit$BUGSoutput$sims.list$deviance[,1]
if(sum(cctfit$BUGSoutput$sims.list$deviance==Inf) >0){
cctfit$BUGSoutput$mean$deviance <- mean(cctfit$BUGSoutput$sims.list$deviance[cctfit$BUGSoutput$sims.list$deviance!=Inf,1],na.rm=TRUE) #Dbar, also known as deviance
cctfit$BUGSoutput$pD <- var(cctfit$BUGSoutput$sims.list$deviance[cctfit$BUGSoutput$sims.list$deviance!=Inf,1],na.rm=TRUE)/2 #pD, variance of the deviance, divided by 2
cctfit$BUGSoutput$DIC <- cctfit$BUGSoutput$mean$deviance + cctfit$BUGSoutput$pD
}else{
cctfit$BUGSoutput$mean$deviance <- mean(cctfit$BUGSoutput$sims.list$deviance) #Dbar, also known as deviance
cctfit$BUGSoutput$pD <- var(cctfit$BUGSoutput$sims.list$deviance[,1])/2 #pD, variance of the deviance, divided by 2
cctfit$BUGSoutput$DIC <- cctfit$BUGSoutput$mean$deviance + cctfit$BUGSoutput$pD
}
}
if(cctfit$whmodel=="LTRM"){
if(cctfit$V==1){
cctfit$V <- 1;
cctfit$BUGSoutput$sims.list$Om <- array(1,c(cctfit$BUGSoutput$n.sims,cctfit$n));
if(length(dim(cctfit$BUGSoutput$sims.list$T)) < 3){
cctfit$BUGSoutput$sims.list$T <- array(cctfit$BUGSoutput$sims.list[["T"]][,],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))}
if(length(dim(cctfit$BUGSoutput$sims.list$gam))<3){cctfit$BUGSoutput$sims.list$gam <- array(cctfit$BUGSoutput$sims.list[["gam"]][,],c(cctfit$BUGSoutput$n.sims,cctfit$C-1,1))}
if(cctfit$itemdiff==TRUE){
if(length(dim(cctfit$BUGSoutput$sims.list$lam)) < 3){
cctfit$BUGSoutput$sims.list$lam <- array(cctfit$BUGSoutput$sims.list[["lam"]][,],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))}
}
}
if(cctfit$itemdiff==FALSE){
if(cctfit$V==1){
cctfit$BUGSoutput$sims.list$lam <- array(1,c(cctfit$BUGSoutput$n.sims,cctfit$m,1))
}
if(cctfit$V > 1){
cctfit$BUGSoutput$sims.list$lam <- array(1,c(cctfit$BUGSoutput$n.sims,cctfit$m,cctfit$V))
}
}
storefulllik <- 0
if(storefulllik==1){
cctfit$BUGSoutput$sims.list$tau <- array(-1,c(cctfit$BUGSoutput$n.sims,cctfit$n,cctfit$m))
cctfit$ppdelta <- array(NA, c(cctfit$n,cctfit$C-1,cctfit$BUGSoutput$n.sims))
cctfit$ppdeltafull <- array(NA, c(cctfit$n,cctfit$C+1,cctfit$BUGSoutput$n.sims))
cctfit$ppdeltafull[,1,] <- -1000000; cctfit$ppdeltafull[,cctfit$C+1,] <- 1000000
cctfit$Lik <- array(NA, c(cctfit$n,cctfit$m,cctfit$BUGSoutput$n.sims));
for(samp in 1:cctfit$BUGSoutput$n.sims){
cctfit$BUGSoutput$sims.list[["tau"]][samp,,] <- cctfit$BUGSoutput$sims.list[["E"]][samp,]*t(1/cctfit$BUGSoutput$sims.list[["lam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])
cctfit$ppdeltafull[,2:(cctfit$C),samp] <- cctfit$BUGSoutput$sims.list[["a"]][samp,]*t(cctfit$BUGSoutput$sims.list[["gam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])+cctfit$BUGSoutput$sims.list[["b"]][samp,]
for(i in 1:cctfit$n){
cctfit$Lik[i,,samp] <- pnorm(cctfit$ppdeltafull[i,cctfit$data[i,]+1,samp] ,cctfit$BUGSoutput$sims.list[["T"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,i]],cctfit$BUGSoutput$sims.list[["tau"]][samp,i,]^-.5)-pnorm(cctfit$ppdeltafull[i,cctfit$data[i,],samp] ,cctfit$BUGSoutput$sims.list[["T"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,i]],cctfit$BUGSoutput$sims.list[["tau"]][samp,i,]^-.5)
if(cctfit$mval==TRUE){cctfit$Lik[cbind(datob$thena,samp)] <- 1}
}
}
if(cctfit$C==2){cctfit$Lik[cctfit$Lik==0] <- 0.001}
cctfit$BUGSoutput$sims.list$deviance <- array(-2*apply(log(cctfit$Lik),3,sum),c(cctfit$BUGSoutput$n.sims,1))
}else{
cctfit$LogLik <- array(NA, c(cctfit$BUGSoutput$n.sims));
Liktmp <- array(NA, c(cctfit$n,cctfit$m))
deltatmp <- array(NA, c(cctfit$n,cctfit$C+1))
deltatmp[,1] <- -100000; deltatmp[,cctfit$C+1] <- 1000000
for(samp in 1:cctfit$BUGSoutput$n.sims){
tautmp <- (cctfit$BUGSoutput$sims.list[["E"]][samp,]*t(exp(cctfit$BUGSoutput$sims.list[["lam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])))^-2
deltatmp[,2:(cctfit$C)] <- cctfit$BUGSoutput$sims.list[["a"]][samp,]*t(cctfit$BUGSoutput$sims.list[["gam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])+cctfit$BUGSoutput$sims.list[["b"]][samp,]
for(i in 1:cctfit$n){
Liktmp[i,] <- pnorm(deltatmp[i,cctfit$data[i,]+1],cctfit$BUGSoutput$sims.list[["T"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,i]],tautmp[i,]^-.5)-pnorm(deltatmp[i,cctfit$data[i,]],cctfit$BUGSoutput$sims.list[["T"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,i]],tautmp[i,]^-.5)
if(cctfit$mval==TRUE){cctfit$Lik[datob$thena] <- 1}
}
if(cctfit$C==2){Liktmp[Liktmp==0] <- 0.001}
cctfit$LogLik[samp] <- sum(log(Liktmp))
}
cctfit$BUGSoutput$sims.list$deviance <- array(-2*(cctfit$LogLik),c(cctfit$BUGSoutput$n.sims,1))
}
cctfit$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="deviance")]]] <- array(cctfit$BUGSoutput$sims.list$deviance,c(cctfit$BUGSoutput$n.keep,cctfit$BUGSoutput$n.chains))
cctfit$BUGSoutput$sims.matrix[,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="deviance")]]] <- cctfit$BUGSoutput$sims.list$deviance[,1]
if(sum(cctfit$BUGSoutput$sims.list$deviance==Inf) >0){
cctfit$BUGSoutput$mean$deviance <- mean(cctfit$BUGSoutput$sims.list$deviance[cctfit$BUGSoutput$sims.list$deviance!=Inf,1],na.rm=TRUE) #Dbar, also known as deviance
cctfit$BUGSoutput$pD <- var(cctfit$BUGSoutput$sims.list$deviance[cctfit$BUGSoutput$sims.list$deviance!=Inf,1],na.rm=TRUE)/2 #pD, variance of the deviance, divided by 2
cctfit$BUGSoutput$DIC <- cctfit$BUGSoutput$mean$deviance + cctfit$BUGSoutput$pD
}else{
cctfit$BUGSoutput$mean$deviance <- mean(cctfit$BUGSoutput$sims.list$deviance) #Dbar, also known as deviance
cctfit$BUGSoutput$pD <- var(cctfit$BUGSoutput$sims.list$deviance[,1])/2 #pD, variance of the deviance, divided by 2
cctfit$BUGSoutput$DIC <- cctfit$BUGSoutput$mean$deviance + cctfit$BUGSoutput$pD
}
}
if(cctfit$whmodel=="CRM"){
if(cctfit$V==1){
cctfit$V <- 1;
cctfit$BUGSoutput$sims.list$Om <- array(1,c(cctfit$BUGSoutput$n.sims,cctfit$n));
if(length(dim(cctfit$BUGSoutput$sims.list$T)) < 3){
cctfit$BUGSoutput$sims.list$T <- array(cctfit$BUGSoutput$sims.list[["T"]][,],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))}
if(cctfit$itemdiff==TRUE){
if(length(dim(cctfit$BUGSoutput$sims.list$lam)) < 3){
cctfit$BUGSoutput$sims.list$lam <- array(cctfit$BUGSoutput$sims.list[["lam"]][,],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))}
}
}
if(cctfit$itemdiff==FALSE){
if(cctfit$V==1){
cctfit$BUGSoutput$sims.list$lam <- array(0,c(cctfit$BUGSoutput$n.sims,cctfit$m,1))
}
if(cctfit$V > 1){
cctfit$BUGSoutput$sims.list$lam <- array(0,c(cctfit$BUGSoutput$n.sims,cctfit$m,cctfit$V))
}
}
storefulllik <- 0
if(storefulllik==1){
cctfit$BUGSoutput$sims.list$tau <- array(-1,c(cctfit$BUGSoutput$n.sims,cctfit$n,cctfit$m))
cctfit$Lik <- array(NA, c(cctfit$n,cctfit$m,cctfit$BUGSoutput$n.sims));
for(samp in 1:cctfit$BUGSoutput$n.sims){
cctfit$BUGSoutput$sims.list[["tau"]][samp,,] <- (cctfit$BUGSoutput$sims.list[["E"]][samp,]*t(exp(cctfit$BUGSoutput$sims.list[["lam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])))^-2
cctfit$Lik[,,samp] <- dnorm(cctfit$data,mean=(cctfit$BUGSoutput$sims.list[["a"]][samp,]*t(cctfit$BUGSoutput$sims.list[["T"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]]))+array(cctfit$BUGSoutput$sims.list[["b"]][samp,],c(cctfit$n,cctfit$m)),sd=cctfit$BUGSoutput$sims.list[["a"]][samp,]*(cctfit$BUGSoutput$sims.list[["tau"]][samp,,]^-.5))
if(cctfit$mval==TRUE){cctfit$Lik[cbind(datob$thena,samp)] <- 1}
}
cctfit$BUGSoutput$sims.list$deviance <- array(-2*apply(log(cctfit$Lik),3,sum),c(cctfit$BUGSoutput$n.sims,1))
}else{
cctfit$LogLik <- array(NA, c(cctfit$BUGSoutput$n.sims));
for(samp in 1:cctfit$BUGSoutput$n.sims){
tautmp <- (cctfit$BUGSoutput$sims.list[["E"]][samp,]*t(exp(cctfit$BUGSoutput$sims.list[["lam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])))^-2
Liktmp <- dnorm(cctfit$data,mean=(cctfit$BUGSoutput$sims.list[["a"]][samp,]*t(cctfit$BUGSoutput$sims.list[["T"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]]))+array(cctfit$BUGSoutput$sims.list[["b"]][samp,],c(cctfit$n,cctfit$m)),sd=cctfit$BUGSoutput$sims.list[["a"]][samp,]*(tautmp^-.5))
if(cctfit$mval==TRUE){Liktmp[datob$thena] <- 1}
cctfit$LogLik[samp] <- sum(log(Liktmp))
}
cctfit$BUGSoutput$sims.list$deviance <- array(-2*(cctfit$LogLik),c(cctfit$BUGSoutput$n.sims,1))
}
cctfit$BUGSoutput$sims.array[,,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="deviance")]]] <- array(cctfit$BUGSoutput$sims.list$deviance,c(cctfit$BUGSoutput$n.keep,cctfit$BUGSoutput$n.chains))
cctfit$BUGSoutput$sims.matrix[,cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="deviance")]]] <- cctfit$BUGSoutput$sims.list$deviance[,1]
if(sum(cctfit$BUGSoutput$sims.list$deviance==Inf) >0){
cctfit$BUGSoutput$mean$deviance <- mean(cctfit$BUGSoutput$sims.list$deviance[cctfit$BUGSoutput$sims.list$deviance!=Inf,1],na.rm=TRUE) #Dbar, also known as deviance
cctfit$BUGSoutput$pD <- var(cctfit$BUGSoutput$sims.list$deviance[cctfit$BUGSoutput$sims.list$deviance!=Inf,1],na.rm=TRUE)/2 #pD, variance of the deviance, divided by 2
cctfit$BUGSoutput$DIC <- cctfit$BUGSoutput$mean$deviance + cctfit$BUGSoutput$pD
}else{
cctfit$BUGSoutput$mean$deviance <- mean(cctfit$BUGSoutput$sims.list$deviance) #Dbar, also known as deviance
cctfit$BUGSoutput$pD <- var(cctfit$BUGSoutput$sims.list$deviance[,1])/2 #pD, variance of the deviance, divided by 2
cctfit$BUGSoutput$DIC <- cctfit$BUGSoutput$mean$deviance + cctfit$BUGSoutput$pD
}
}
message(paste("DIC : ",round(cctfit$BUGSoutput$DIC,2)," pD : ",round(cctfit$BUGSoutput$pD,2),sep=""))
if(gui==TRUE){
guidat <- get("guidat", pkg_globals)
tkconfigure(guidat$plotresults.but, state="normal")
tkconfigure(guidat$doppc.but, state="normal")
tkconfigure(guidat$exportresults.but, state="normal")
tkconfigure(guidat$printfit.but, state="normal")
tkconfigure(guidat$summary.but, state="normal")
tkconfigure(guidat$sendconsole.but, state="normal")
tkconfigure(guidat$memb.but, state="normal")
tkconfigure(guidat$mvest.but, state="normal")
#if(guidat$mval == TRUE){tkconfigure(guidat$mvest.but, state="normal")}
tkconfigure(guidat$traceplotdiscrete.but, state="normal")
tkconfigure(guidat$traceplotall.but, state="normal")
cctfit$guidat <- guidat
assign("guidat", guidat, pkg_globals)
class(cctfit) <- c("cct",class(cctfit))
return(cctfit)
}else{
class(cctfit) <- c("cct",class(cctfit))
return(cctfit)
}
}
######################
#Function for the 'Plot Results' Button
#- Creates a sophisticated plot of the posterior results, which includes
# the posterior means of each parameters and their highest density intervals (HDIs, see Kruschke 2011)
#- Denotes a different symbol for the parameters pertaining to each culture
#- This function is also called during the file export, and .eps and .jpeg's of the plot are saved
######################
plotresultsfuncbutton <- function() {
guidat <- get("guidat", pkg_globals);
plotresultsfunc(cctfit=get(guidat$varname, pkg_globals),gui=TRUE)
}
plotresultsfunc <- function(cctfit,saveplots=0,savedir="",gui=FALSE) {
hdi <- function(sampleVec, credMass=0.95){
sortedPts = sort(sampleVec)
ciIdxInc = floor(credMass * length( sortedPts) )
nCIs = length( sortedPts) - ciIdxInc
ciwidth = rep(0, nCIs)
for(i in 1:nCIs){
ciwidth[i] = sortedPts[i + ciIdxInc] - sortedPts[i]}
HDImin = sortedPts[which.min(ciwidth)]
HDImax = sortedPts[which.min(ciwidth)+ciIdxInc]
HDIlim = c(HDImin,HDImax)
return(HDIlim)
}
Mode <- function(x) {ux <- unique(x); ux[which.max(tabulate(match(x, ux)))] }
cctfit$respmem <- apply(cctfit$BUGSoutput$sims.list$Om[,],2,Mode)
if(saveplots==1){jpeg(file.path(gsub(".Rdata","results.jpg",savedir)),width = 6, height = 6, units = "in", pointsize = 12,quality=100,res=400)}
if(saveplots==2){postscript(file=file.path(gsub(".Rdata","results.eps",savedir)), onefile=FALSE, horizontal=FALSE, width = 6, height = 6, paper="special", family="Times")}
if(cctfit$whmodel=="GCM"){
par(oma=c(0,0,0,0),mar=c(4,4,3,1),mgp=c(2.25,.75,0),mfrow=c(2,2))
sym <- c(21,22, 23, 24, 25, 3, 4, 5, 7, 8, 9, 10, 11, 12, 13, 14) # possible symbols:
bgcol <- c("black",NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
#if(cctfit$V==1){bgcol <- c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)}
if(cctfit$V==1){
if(length(dim(cctfit$BUGSoutput$sims.list[["Z"]])) < 3){
cctfit$BUGSoutput$sims.list[["Z"]] <- array(cctfit$BUGSoutput$sims.list[["Z"]],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))}
plot(1:cctfit$m,rep(NA,cctfit$m),main=expression(paste("Item Truth (",Z[vk],")")),xlab="Item",ylab="Posterior Mean Value",las=1,ylim=c(0,1),pch=21,bg="white")
for(i in 1:dim(cctfit$BUGSoutput$sims.list[["Z"]])[3]){
points(apply(cctfit$BUGSoutput$sims.list[["Z"]][,,i],c(2),mean),xlab=expression(paste("Item Truth (",Z[vk],") Per Cluster")),ylab="Posterior Mean Value",las=1,ylim=c(0,1),pch=sym[i],bg=bgcol[i])
}}else{
plot(1:cctfit$m,rep(NA,cctfit$m),main=expression(paste("Item Truth (",Z[vk],") Per Culture")),xlab="Item",ylab="Posterior Mean Value",las=1,ylim=c(0,1),pch=21,bg="white")
for(i in 1:dim(cctfit$BUGSoutput$sims.list[["Z"]])[3]){
points(apply(cctfit$BUGSoutput$sims.list[["Z"]][,,i],c(2),mean),xlab=expression(paste("Item Truth (",Z[vk],") Per Cluster")),ylab="Posterior Mean Value",las=1,ylim=c(0,1),pch=sym[i],bg=bgcol[i])
}
}
if(cctfit$itemdiff==TRUE){
if(cctfit$V==1){
if(length(dim(cctfit$BUGSoutput$sims.list[["lam"]])) < 3){
cctfit$BUGSoutput$sims.list[["lam"]] <- array(cctfit$BUGSoutput$sims.list[["lam"]],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))}
hditmp <- apply(cctfit$BUGSoutput$sims.list[["lam"]][,,1],2,hdi)
plot(1:cctfit$m,rep(NA,cctfit$m),main=expression(paste("Item Difficulty (",lambda[vk],")")),xlab="Item",ylim=c(0,1),ylab="Posterior Mean Value",las=1,pch=21,bg="white")
for(i in 1:dim(cctfit$BUGSoutput$sims.list[["lam"]])[3]){
points(apply(cctfit$BUGSoutput$sims.list[["lam"]][,,i],c(2),mean),xlab=expression(paste("Item Difficulty (",lambda[vk],") Per Cluster")),ylab="Posterior Mean Value",las=1,ylim=c(0,1),pch=sym[i],bg=bgcol[i])
}
segments(1:cctfit$m,hditmp[1,], 1:cctfit$m, hditmp[2,])
arrows(1:cctfit$m,hditmp[1,], 1:cctfit$m, hditmp[2,],code=3,angle=90,length=.025)
}else{plot(1:cctfit$m,rep(NA,cctfit$m),main=expression(paste("Item Difficulty (",lambda[vk],") Per Culture")),xlab="Item",ylim=c(0,1),ylab="Posterior Mean Value",las=1,pch=21)
for(i in 1:dim(cctfit$BUGSoutput$sims.list[["lam"]])[3]){
points(apply(cctfit$BUGSoutput$sims.list[["lam"]][,,i],c(2),mean),xlab=expression(paste("Item Difficulty (",lambda[vk],") Per Cluster")),ylab="Posterior Mean Value",las=1,ylim=c(0,1),pch=sym[i],bg=bgcol[i])
}
}
}else{
plot(c(1:cctfit$m),rep(.5,cctfit$m),main=expression(paste("Item Difficulty (",lambda[vk],")")),xlab="Item",ylab="Posterior Mean Value",las=1,ylim=c(0,1),pch=sym[1],bg="white",col="grey")
points(c(par("usr")[1],par("usr")[2]),c(par("usr")[3],par("usr")[4]),type="l",col="grey")
points(c(par("usr")[1],par("usr")[2]),c(par("usr")[4],par("usr")[3]),type="l",col="grey")
}
plot(cctfit$BUGSoutput$mean$th,main=expression(paste("Respondent Competency (",theta[i],")")),xlab="Respondent",ylab="Posterior Mean Value",las=1,ylim=c(0,1),pch=sym[cctfit$respmem],bg=bgcol[cctfit$respmem])
hditmp <- apply(cctfit$BUGSoutput$sims.list$th,2,hdi)
segments(1:cctfit$n,hditmp[1,], 1:cctfit$n, hditmp[2,])
arrows(1:cctfit$n,hditmp[1,], 1:cctfit$n, hditmp[2,],code=3,angle=90,length=.025)
plot(cctfit$BUGSoutput$mean$g,main=expression(paste("Respondent Guessing Bias (",g[i],")")),xlab="Respondent",ylab="Posterior Mean Value",las=1,ylim=c(0,1),pch=sym[cctfit$respmem],bg=bgcol[cctfit$respmem])
hditmp <- apply(cctfit$BUGSoutput$sims.list$g,2,hdi)
segments(1:cctfit$n,hditmp[1,], 1:cctfit$n, hditmp[2,])
arrows(1:cctfit$n,hditmp[1,], 1:cctfit$n, hditmp[2,],code=3,angle=90,length=.025)
}
if(cctfit$whmodel=="LTRM"){
invlogit <- function(x){x <- 1 / (1 + exp(-x)); return(x)}
par(oma=c(0,0,0,0),mar=c(4,4,3,1),mgp=c(2.25,.75,0),mfrow=c(2,2))
sym <- c(21,22, 23, 24, 25, 3, 4, 5, 7, 8, 9, 10, 11, 12, 13, 14) # possible symbols:
bgcol <- c("black",NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
if(cctfit$C==2){useinvlogit <- TRUE}else{useinvlogit <- FALSE}
if(useinvlogit==TRUE){
#tmp1 <- cctfit$BUGSoutput$sims.list[["T"]]; cctfit$BUGSoutput$sims.list[["T"]] <- invlogit(cctfit$BUGSoutput$sims.list[["T"]])
#tmp2 <- cctfit$BUGSoutput$sims.list[["gam"]]; cctfit$BUGSoutput$sims.list[["gam"]] <- invlogit(cctfit$BUGSoutput$sims.list[["gam"]])
#tmp3 <- cctfit$BUGSoutput$sims.list[["b"]]; cctfit$BUGSoutput$sims.list[["b"]] <- invlogit(cctfit$BUGSoutput$sims.list[["b"]])
tmp1 <- cctfit$BUGSoutput$sims.list[["T"]]; cctfit$BUGSoutput$sims.list[["T"]] <- invlogit(cctfit$BUGSoutput$sims.list[["T"]])
tmp2 <- cctfit$BUGSoutput$sims.list[["gam"]]; cctfit$BUGSoutput$sims.list[["gam"]] <- invlogit(cctfit$BUGSoutput$sims.list[["gam"]])
tmp3 <- cctfit$BUGSoutput$sims.list[["b"]]; cctfit$BUGSoutput$sims.list[["b"]] <- invlogit(cctfit$BUGSoutput$sims.list[["b"]])
}
if(cctfit$V==1){
newm <- ceiling(cctfit$m*1.092)
if(length(dim(cctfit$BUGSoutput$sims.list[["T"]])) <3){
cctfit$BUGSoutput$sims.list[["T"]] <- array(cctfit$BUGSoutput$sims.list[["T"]],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))
}
if(length(dim(cctfit$BUGSoutput$sims.list[["gam"]])) <3){
cctfit$BUGSoutput$sims.list[["gam"]] <- array(cctfit$BUGSoutput$sims.list[["gam"]],c(cctfit$BUGSoutput$n.sims,cctfit$C-1,1))
}
hditmp <- apply(cctfit$BUGSoutput$sims.list[["T"]][,,1],2,hdi)
if(useinvlogit==TRUE){
plot(1:newm,rep(NA,newm),main=expression(paste("Item Truth (",T[vk],") and Thresholds (",gamma[vc],")")),xlab="Item",ylim=c(0,1),ylab="Posterior Mean Invlogit Value",las=1,pch=21,bg="white",axes=FALSE)
}else{plot(1:newm,rep(NA,newm),main=expression(paste("Item Truth (",T[vk],") and Thresholds (",gamma[vc],")")),xlab="Item",ylim=c(min(hditmp,min(apply(cctfit$BUGSoutput$sims.list[["gam"]],c(2,3),mean))),max(hditmp,max(apply(cctfit$BUGSoutput$sims.list[["gam"]],c(2,3),mean)))),ylab="Posterior Mean Value",las=1,pch=21,bg="white",axes=FALSE)
}
for(i in 1:dim(cctfit$BUGSoutput$sims.list[["T"]])[3]){
points(apply(cctfit$BUGSoutput$sims.list[["T"]][,,i],c(2),mean),xlab=expression(paste("Item Truth (",T[vk],") Per Cluster")),ylab="Posterior Mean Value",las=1,ylim=c(0,1),pch=sym[i],bg=bgcol[i])
if(cctfit$C==2){
text(newm,mean(cctfit$BUGSoutput$sims.list[["gam"]][,,i]),labels=sapply(c(1:(cctfit$C-1)),function(x) as.expression(substitute(list(gamma[x]),list(x=x)))))
segments(1,mean(cctfit$BUGSoutput$sims.list[["gam"]][,,i]),max(cctfit$m+1,newm-2),mean(cctfit$BUGSoutput$sims.list[["gam"]][,,i]),lty=2);
}
else{
text(newm,apply(cctfit$BUGSoutput$sims.list[["gam"]][,,i],c(2),mean),labels=sapply(c(1:(cctfit$C-1)),function(x) as.expression(substitute(list(gamma[x]),list(x=x)))))
segments(1,apply(cctfit$BUGSoutput$sims.list[["gam"]][,,i],c(2),mean),max(cctfit$m+1,newm-2),apply(cctfit$BUGSoutput$sims.list[["gam"]][,,i],c(2),mean),lty=2);
}
segments(1:cctfit$m,hditmp[1,], 1:cctfit$m, hditmp[2,])
arrows(1:cctfit$m,hditmp[1,], 1:cctfit$m, hditmp[2,],code=3,angle=90,length=.025)
box()
axis(2, labels = TRUE,las=1)
axis(side = 1,labels=TRUE)
axis(side = 1, at = newm, labels = expression(gamma[c]) )
}}else{
newm <- ceiling(cctfit$m*1.14)
if(cctfit$C==2){
if(useinvlogit==TRUE){
plot(1:newm,rep(NA,newm),main=expression(paste("Item Truth (",T[vk],") and Thresholds (",gamma[vc],") Per Culture")),xlab="Item",ylim=c(0,1),ylab="Posterior Mean Value",las=1,pch=21,bg="white",axes=FALSE)
}else{plot(1:newm,rep(NA,newm),main=expression(paste("Item Truth (",T[vk],") and Thresholds (",gamma[vc],") Per Culture")),xlab="Item",ylim=c(min(min(apply(cctfit$BUGSoutput$sims.list[["T"]],c(2,3),mean))
,min(apply(cctfit$BUGSoutput$sims.list[["gam"]],c(3),mean))),max(max(apply(cctfit$BUGSoutput$sims.list[["T"]],c(2,3),mean))
,max(apply(cctfit$BUGSoutput$sims.list[["gam"]],c(3),mean)))),ylab="Posterior Mean Value",las=1,pch=21,bg="white",axes=FALSE)
}
for(i in 1:dim(cctfit$BUGSoutput$sims.list[["T"]])[3]){
points(apply(cctfit$BUGSoutput$sims.list[["T"]][,,i],c(2),mean),xlab=expression(paste("Item Truth (",T[vk],") Per Cluster")),ylab="Posterior Mean Value",las=1,ylim=c(0,1),pch=sym[i],bg=bgcol[i])
text(cctfit$m+(((newm-cctfit$m)/dim(cctfit$BUGSoutput$sims.list[["T"]])[3])*i),mean(cctfit$BUGSoutput$sims.list[["gam"]][,,i]),labels=sapply(c(1:(cctfit$C-1)),function(x) as.expression(substitute(list(gamma[x]),list(x=x)))))
}
}else{
if(useinvlogit==TRUE){
plot(1:newm,rep(NA,newm),main=expression(paste("Item Truth (",T[vk],") and Thresholds (",gamma[vc],") Per Culture")),xlab="Item",ylim=c(0,1),ylab="Posterior Mean Value",las=1,pch=21,bg="white",axes=FALSE)
}else{
plot(1:newm,rep(NA,newm),main=expression(paste("Item Truth (",T[vk],") and Thresholds (",gamma[vc],") Per Culture")),xlab="Item",ylim=c(min(min(apply(cctfit$BUGSoutput$sims.list[["T"]],c(2,3),mean))
,min(apply(cctfit$BUGSoutput$sims.list[["gam"]][,,],c(2,3),mean))),max(max(apply(cctfit$BUGSoutput$sims.list[["T"]],c(2,3),mean))
,max(apply(cctfit$BUGSoutput$sims.list[["gam"]][,,],c(2,3),mean)))),ylab="Posterior Mean Value",las=1,pch=21,bg="white",axes=FALSE)
}
for(i in 1:dim(cctfit$BUGSoutput$sims.list[["T"]])[3]){
points(apply(cctfit$BUGSoutput$sims.list[["T"]][,,i],c(2),mean),xlab=expression(paste("Item Truth (",T[vk],") Per Cluster")),ylab="Posterior Mean Value",las=1,ylim=c(0,1),pch=sym[i],bg=bgcol[i])
text(cctfit$m+(((newm-cctfit$m)/dim(cctfit$BUGSoutput$sims.list[["T"]])[3])*i),apply(cctfit$BUGSoutput$sims.list[["gam"]][,,i],c(2),mean),labels=sapply(c(1:(cctfit$C-1)),function(x) as.expression(substitute(list(gamma[x]),list(x=x)))))
}
}
box()
axis(2, labels = TRUE,las=1)
axis(side=1,at=axTicks(side=1)[axTicks(side=1)<= cctfit$m])
axis(side = 1, at = cctfit$m+(((newm-cctfit$m)/dim(cctfit$BUGSoutput$sims.list[["T"]])[3])*(1:2)),
labels = sapply(1:2,function(x) as.expression(substitute(list(gamma[x*c]),list(x=x)))) )
}
if(cctfit$itemdiff==TRUE){
if(cctfit$V==1){
if(length(dim(cctfit$BUGSoutput$sims.list[["lam"]])) < 3){
cctfit$BUGSoutput$sims.list[["lam"]] <- array(cctfit$BUGSoutput$sims.list[["lam"]],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))}
hditmp <- apply(cctfit$BUGSoutput$sims.list[["lam"]][,,1],2,hdi)
plot(1:cctfit$m,rep(NA,cctfit$m),main=expression(paste("Item Difficulty (",lambda[vk],")")),xlab="Item",ylim=c(min(hditmp),max(hditmp)),ylab="Posterior Mean Log Value",las=1,pch=21,bg="white")
for(i in 1:dim(cctfit$BUGSoutput$sims.list[["lam"]])[3]){
points(apply(cctfit$BUGSoutput$sims.list[["lam"]][,,i],c(2),mean),xlab=expression(paste("Item Difficulty (",lambda[vk],") Per Cluster")),ylab="Posterior Mean Log Value",las=1,ylim=c(0,1),pch=sym[i],bg=bgcol[i])
}
segments(1:cctfit$m,hditmp[1,], 1:cctfit$m, hditmp[2,])
arrows(1:cctfit$m,hditmp[1,], 1:cctfit$m, hditmp[2,],code=3,angle=90,length=.025)
}else{plot(1:cctfit$m,rep(NA,cctfit$m),main=expression(paste("Item Difficulty (",lambda[vk],") Per Culture")),xlab="Item",ylim=c(min(apply(cctfit$BUGSoutput$sims.list[["lam"]],c(2,3),mean)),max(apply(cctfit$BUGSoutput$sims.list[["lam"]],c(2,3),mean))),ylab="Posterior Mean Value",las=1,pch=21,bg="white")
for(i in 1:dim(cctfit$BUGSoutput$sims.list[["lam"]])[3]){
points(apply(cctfit$BUGSoutput$sims.list[["lam"]][,,i],c(2),mean),xlab=expression(paste("Item Difficulty (",lambda[vk],") Per Cluster")),ylab="Posterior Mean Log Value",las=1,ylim=c(0,1),pch=sym[i],bg=bgcol[i])
}
}
}else{
plot(c(1:cctfit$m),rep(.5,cctfit$m),main=expression(paste("Item Difficulty (",lambda[vk],")")),xlab="Item",ylab="Posterior Mean Log Value",las=1,ylim=c(0,1),pch=sym[1],bg="white",col="grey")
points(c(par("usr")[1],par("usr")[2]),c(par("usr")[3],par("usr")[4]),type="l",col="grey")
points(c(par("usr")[1],par("usr")[2]),c(par("usr")[4],par("usr")[3]),type="l",col="grey")
}
hditmp <- apply(cctfit$BUGSoutput$sims.list$E,2,hdi)
plot(cctfit$BUGSoutput$mean$E,main=expression(paste("Respondent Standard Error (",E[i],")")),xlab="Respondent",ylab="Posterior Mean Value",las=1,ylim=c(0,max(hditmp)),pch=sym[cctfit$respmem],bg=bgcol[cctfit$respmem])
segments(1:cctfit$n,hditmp[1,], 1:cctfit$n, hditmp[2,])
arrows(1:cctfit$n,hditmp[1,], 1:cctfit$n, hditmp[2,],code=3,angle=90,length=.025)
if(useinvlogit==TRUE){
plot(invlogit(cctfit$BUGSoutput$mean$b)-.5,main=expression(paste("Respondent Shift Bias (",b[i],")")),xlab="Respondent",ylab="Posterior Mean Value",las=1,ylim=c(-.5,.5),pch=sym[cctfit$respmem],bg=bgcol[cctfit$respmem])
hditmp <- apply(cctfit$BUGSoutput$sims.list$b-.5,2,hdi) # because this is invprobit transformed still
segments(1:cctfit$n,hditmp[1,], 1:cctfit$n, hditmp[2,])
arrows(1:cctfit$n,hditmp[1,], 1:cctfit$n, hditmp[2,],code=3,angle=90,length=.025)
}else{
loga <- apply(log(cctfit$BUGSoutput$sims.list[["a"]]),2,mean)
plot(-5000, -5000, main=expression(paste("Category Usage Bias (",a[i]," and ",b[i],")")),xlim=c(min(-1.6,-abs(cctfit$BUGSoutput$mean$b)),max(1.6,cctfit$BUGSoutput$mean$b)), ylim=c(min(-1,loga),max(1,loga)),xlab=expression(paste("Respondent Shift Bias (",b[i],")")),
ylab=expression(paste("Log Respondent Scale Bias (",a[i],")")),las=1,pch=21)
points(cbind(cctfit$BUGSoutput$mean$b,loga),xlab="Respondent",ylab="Posterior Mean Value",las=1,ylim=c(0,max(hditmp)),pch=sym[cctfit$respmem],bg=bgcol[cctfit$respmem])
segments(.625*par("usr")[1],0,.625*par("usr")[2],0)
segments(0,.725*par("usr")[3],0,.725*par("usr")[4])
text(0,.875*par("usr")[3],"Middle Categories",cex=.8)
text(0,.875*par("usr")[4],"Outer Categories",cex=.8)
text(.775*par("usr")[1],0,"Left \n Categories",cex=.8)
text(.775*par("usr")[2],0,"Right \n Categories",cex=.8)
}
if(useinvlogit==TRUE){
cctfit$BUGSoutput$sims.list[["T"]] <- tmp1
cctfit$BUGSoutput$sims.list[["gam"]] <- tmp2
cctfit$BUGSoutput$sims.list[["b"]] <- tmp3
}
}
if(cctfit$whmodel=="CRM"){
par(oma=c(0,0,0,0),mar=c(4,4,3,1),mgp=c(2.25,.75,0),mfrow=c(2,2))
sym <- c(21,22, 23, 24, 25, 3, 4, 5, 7, 8, 9, 10, 11, 12, 13, 14) # possible symbols:
bgcol <- c("black",NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
if(cctfit$V==1){
if(length(dim(cctfit$BUGSoutput$sims.list[["T"]])) < 3){
cctfit$BUGSoutput$sims.list[["T"]] <- array(cctfit$BUGSoutput$sims.list[["T"]],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))}
hditmp <- apply(cctfit$BUGSoutput$sims.list[["T"]][,,1],2,hdi)
plot(1:cctfit$m,rep(NA,cctfit$m),main=expression(paste("Item Truth (",T[vk],")")),xlab="Item",ylim=c(min(hditmp),max(hditmp)),ylab="Posterior Mean Value",las=1,pch=21,bg="white")
}else{plot(1:cctfit$m,rep(NA,cctfit$m),main=expression(paste("Item Truth (",T[vk],") Per Culture")),xlab="Item",ylim=c(min(apply(cctfit$BUGSoutput$sims.list[["T"]],c(2,3),mean)),max(apply(cctfit$BUGSoutput$sims.list[["T"]],c(2,3),mean))),ylab="Posterior Mean Value",las=1,pch=21,bg="white")}
for(i in 1:dim(cctfit$BUGSoutput$sims.list[["T"]])[3]){
points(apply(cctfit$BUGSoutput$sims.list[["T"]][,,i],c(2),mean),xlab=expression(paste("Item Truth (",T[vk],") Per Cluster")),ylab="Posterior Mean Value",las=1,ylim=c(0,1),pch=sym[i],bg=bgcol[i])
}
if(cctfit$V==1){
segments(1:cctfit$m,hditmp[1,], 1:cctfit$m, hditmp[2,])
arrows(1:cctfit$m,hditmp[1,], 1:cctfit$m, hditmp[2,],code=3,angle=90,length=.025)
}
if(cctfit$itemdiff==TRUE){
if(cctfit$V==1){
if(length(dim(cctfit$BUGSoutput$sims.list[["lam"]])) < 3){
cctfit$BUGSoutput$sims.list[["lam"]] <- array(cctfit$BUGSoutput$sims.list[["lam"]],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))}
hditmp <- apply(cctfit$BUGSoutput$sims.list[["lam"]][,,1],2,hdi)
plot(1:cctfit$m,rep(NA,cctfit$m),main=expression(paste("Item Difficulty (",lambda[vk],")")),xlab="Item",ylim=c(min(hditmp),max(hditmp)),ylab="Posterior Mean Value",las=1,pch=21,bg="white")
for(i in 1:dim(cctfit$BUGSoutput$sims.list[["lam"]])[3]){
points(apply(cctfit$BUGSoutput$sims.list[["lam"]][,,i],c(2),mean),xlab=expression(paste("Item Difficulty (",lambda[vk],") Per Cluster")),ylab="Posterior Mean Value",las=1,ylim=c(0,1),pch=sym[i],bg=bgcol[i])
}
segments(1:cctfit$m,hditmp[1,], 1:cctfit$m, hditmp[2,])
arrows(1:cctfit$m,hditmp[1,], 1:cctfit$m, hditmp[2,],code=3,angle=90,length=.025)
}else{plot(1:cctfit$m,rep(NA,cctfit$m),main=expression(paste("Item Difficulty (",lambda[vk],") Per Culture")),xlab="Item",ylim=c(min(apply(cctfit$BUGSoutput$sims.list[["lam"]],c(2,3),mean)),max(apply(cctfit$BUGSoutput$sims.list[["lam"]],c(2,3),mean))),ylab="Posterior Mean Log Value",las=1,pch=21,bg="white")
for(i in 1:dim(cctfit$BUGSoutput$sims.list[["lam"]])[3]){
points(apply(cctfit$BUGSoutput$sims.list[["lam"]][,,i],c(2),mean),xlab=expression(paste("Item Difficulty (",lambda[vk],") Per Cluster")),ylab="Posterior Mean Log Value",las=1,ylim=c(0,1),pch=sym[i],bg=bgcol[i])
}
}
}else{
plot(c(1:cctfit$m),rep(.5,cctfit$m),main=expression(paste("Item Difficulty (",lambda[vk],")")),xlab="Item",ylab="Posterior Mean Log Value",las=1,ylim=c(0,1),pch=sym[1],bg="white",col="grey")
points(c(par("usr")[1],par("usr")[2]),c(par("usr")[3],par("usr")[4]),type="l",col="grey")
points(c(par("usr")[1],par("usr")[2]),c(par("usr")[4],par("usr")[3]),type="l",col="grey")
}
hditmp <- apply(cctfit$BUGSoutput$sims.list$E,2,hdi)
plot(cctfit$BUGSoutput$mean$E,main=expression(paste("Respondent Standard Error (",E[i],")")),xlab="Respondent",ylab="Posterior Mean Value",las=1,ylim=c(min(hditmp),max(hditmp)),pch=sym[cctfit$respmem],bg=bgcol[cctfit$respmem])
segments(1:cctfit$n,hditmp[1,], 1:cctfit$n, hditmp[2,])
arrows(1:cctfit$n,hditmp[1,], 1:cctfit$n, hditmp[2,],code=3,angle=90,length=.025)
#text(cctfit$n/4.25,1.3*par("usr")[3],"More Knowledgeable",cex=.8)
#text(cctfit$n/4.25,.9*par("usr")[4],"Less Knowledgeable",cex=.8)
loga <- apply(log(cctfit$BUGSoutput$sims.list[["a"]]),2,mean)
plot(-5000, -5000, main=expression(paste("Category Usage Bias (",a[i]," and ",b[i],")")),xlim=c(min(-1.6,-abs(cctfit$BUGSoutput$mean$b)),max(1.6,cctfit$BUGSoutput$mean$b)), ylim=c(min(-1,loga),max(1,loga)),xlab=expression(paste("Respondent Shift Bias (",b[i],")")),
ylab=expression(paste("Log Respondent Scale Bias (",a[i],")")),las=1,pch=21)
points(cbind(cctfit$BUGSoutput$mean$b,loga),xlab="Respondent",ylab="Posterior Mean Value",las=1,ylim=c(0,max(hditmp)),pch=sym[cctfit$respmem],bg=bgcol[cctfit$respmem])
segments(.625*par("usr")[1],0,.625*par("usr")[2],0)
segments(0,.725*par("usr")[3],0,.725*par("usr")[4])
text(0,.875*par("usr")[3],"Middle Categories",cex=.8)
text(0,.875*par("usr")[4],"Outer Categories",cex=.8)
text(.775*par("usr")[1],0,"Left \n Categories",cex=.8)
text(.775*par("usr")[2],0,"Right \n Categories",cex=.8)
}
if(saveplots==1 || saveplots==2){dev.off()}
saveplots <- 0
}
######################
#Function for the 'Run Checks' Button
#- Calculates the posterior predictive data for the model that was applied
#- Calculates 2 important posterior predictive checks and plots them
#- The VDI check percentiles are reported in the R console, and saved to cctfit$VDIperc
#- Performs factor analysis using fa() and polychoric() (if selected)
#- This function is also called in the file export and saves .eps and .jpegs of the plot
######################
ppcfuncbutton <- function(){
guidat <- get("guidat", pkg_globals)
cctfit <- get(guidat$varname, pkg_globals)
recalc <- FALSE
if(cctfit$whmodel=="LTRM" && cctfit$checksrun==TRUE){
if(cctfit$polycor!=as.logical(as.numeric(tclvalue(guidat$polyvar)))){
recalc <- TRUE
}
}
assign(guidat$varname,
ppcfunc(cctfit=cctfit,gui=TRUE,rerunchecks=recalc,polych=as.logical(as.numeric(tclvalue(guidat$polyvar)))),
pkg_globals
)
}
ppcfunc <- function(cctfit,saveplots=0,savedir="",gui=FALSE,polych=FALSE,rerunchecks=FALSE,doplot=TRUE) {
if(gui==TRUE){guidat <- get("guidat", pkg_globals)}
if(cctfit$checksrun==FALSE || rerunchecks==TRUE){
message("\n ...One moment, calculating posterior predictive checks")
usesubset <- 1
if(usesubset==1){
subsetsize <- min(500,cctfit$BUGSoutput$n.sims)
indices <- sample(cctfit$BUGSoutput$n.sims,min(subsetsize,cctfit$BUGSoutput$n.sims))
if(cctfit$whmodel=="GCM"){
if(cctfit$V==1){
cctfit$V <- 1;
cctfit$BUGSoutput$sims.list$Om <- array(1,c(cctfit$BUGSoutput$n.sims,cctfit$n));
if(length(dim(cctfit$BUGSoutput$sims.list$Z)) < 3){
cctfit$BUGSoutput$sims.list$Z <- array(cctfit$BUGSoutput$sims.list[["Z"]][,],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))}
if(cctfit$itemdiff==TRUE){
if(length(dim(cctfit$BUGSoutput$sims.list$lam)) < 3){
cctfit$BUGSoutput$sims.list$lam <- array(cctfit$BUGSoutput$sims.list[["lam"]][,],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))}
}
}
if(cctfit$itemdiff==FALSE){
if(cctfit$V==1){
cctfit$BUGSoutput$sims.list$lam <- array(.5,c(cctfit$BUGSoutput$n.sims,cctfit$m,1))
}
if(cctfit$V > 1){
cctfit$BUGSoutput$sims.list$lam <- array(.5,c(cctfit$BUGSoutput$n.sims,cctfit$m,cctfit$V))
} }
cctfit$ppY <- array(NA, c(cctfit$n,cctfit$m,subsetsize));
for(samp in indices){
tautmp <- (cctfit$BUGSoutput$sims.list[["th"]][samp,]*t(1-cctfit$BUGSoutput$sims.list[["lam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])) /
( (cctfit$BUGSoutput$sims.list[["th"]][samp,]*t(1-cctfit$BUGSoutput$sims.list[["lam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])) +
((1-cctfit$BUGSoutput$sims.list[["th"]][samp,])*t(cctfit$BUGSoutput$sims.list[["lam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])) )
cctfit$ppY[,,which(indices==samp)] <- matrix(rbinom((cctfit$n*cctfit$m),1,
(t(cctfit$BUGSoutput$sims.list[["Z"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])*tautmp ) +
t(t(1-tautmp)%*%diag(cctfit$BUGSoutput$sims.list[["g"]][samp,])) ), cctfit$n,cctfit$m)
}
if(cctfit$mval==TRUE){
for(i in 1:dim(cctfit$datob$thena)[1]){cctfit$data[cctfit$datob$thena[i,1],cctfit$datob$thena[i,2]] <- Mode(cctfit$ppY[cctfit$datob$thena[i,1],cctfit$datob$thena[i,2],])}
cctfit$MVest <- cbind(cctfit$datob$thena,cctfit$data[cctfit$datob$thena])
colnames(cctfit$MVest) <- c("Pers","Item","Resp")
}
cctfit$mean$ppY <- rowMeans(cctfit$ppY[,,],dims=2)
leigv <- 12
options(warn=-3)
ind <- indices
eigv <- matrix(-1,length(ind),leigv)
tmp <- apply(cctfit$ppY,c(1,3),function(x) t(x))
tmp2 <- apply(tmp[,,],3,function(x) cor(x))
tmp3 <- array(tmp2,c(cctfit$n,cctfit$n,subsetsize))
for(i in 1:length(ind)){
suppressMessages(try(eigv[i,] <- fa(tmp3[,,i])$values[1:leigv],silent=TRUE))}
wch <- -which(eigv[,1]==-1);
if(length(wch)==0){cctfit$ppeig <- eigv}else{cctfit$ppeig <- eigv[-which(eigv[,1]==-1),]}
if(sum(apply(cctfit$data,1,function(x) sd(x,na.rm=TRUE))==0) > 0){
tmp <- cctfit$data
if(sum(apply(cctfit$datob$dat,1,function(x) sd(x,na.rm=TRUE))==0)==1){
tmp[which(apply(tmp,1,function(x) sd(x,na.rm=TRUE))==0),1] <- min(tmp[which(apply(tmp,1,function(x) sd(x,na.rm=TRUE))==0),1]+.01,.99)
}else{
tmp[which(apply(tmp,1,function(x) sd(x,na.rm=TRUE))==0),1] <- sapply(apply(tmp[which(apply(tmp,1,function(x) sd(x,na.rm=TRUE))==0),],1,mean),function(x) min(x+.01,.99))
}
cctfit$dateig <- suppressMessages(fa(cor(t(tmp)))$values[1:leigv])
}else{
cctfit$dateig <- suppressMessages(fa(cor(t(cctfit$data)))$values[1:leigv])
}
options(warn=0)
if(cctfit$V==1){
varvec <- matrix(-1,length(ind),dim(cctfit$ppY)[2])
vdi <- matrix(-1,dim(cctfit$ppY)[3],1)
for(i in 1:length(ind)){
varvec[i,] <- apply(cctfit$ppY[,,i],2,var)
}
vdi <- apply(varvec,1,var)
cctfit$ppVDI <- vdi
cctfit$datVDI <- var(apply(cctfit$data,2,var))
}
options(warn=-3)
if(cctfit$V>1){
varvec <- array(NA,c(length(ind),dim(cctfit$ppY)[2],cctfit$V))
vdi <- matrix(NA,dim(cctfit$ppY)[3],cctfit$V)
cctfit$datVDI <- array(-1,c(1,cctfit$V))
for(v in 1:cctfit$V){
for(i in 1:length(ind)){
if(sum(cctfit$BUGSoutput$sims.list$Om[ind[i],]==v)>1){
suppressMessages(try(varvec[i,,v] <- apply(cctfit$ppY[cctfit$BUGSoutput$sims.list$Om[ind[i],]==v,,i],2,var),silent=TRUE))
}else{
suppressMessages(try(varvec[i,,v] <- var(cctfit$ppY[cctfit$BUGSoutput$sims.list$Om[ind[i],]==v,,i]),silent=TRUE))
}
}
if(sum(cctfit$respmem==v)>1){
cctfit$datVDI[v] <- var(apply(cctfit$data[cctfit$respmem==v,],2,var))
}else{
cctfit$datVDI[v] <- var(cctfit$data[cctfit$respmem==v,])
}
}
vdi <- apply(varvec,c(1,3),var)
cctfit$ppVDI <- vdi
options(warn=0)
if(sum(apply(cctfit$ppVDI,2,function(x) all(x==0,na.rm=TRUE)))>=1){
for(i in which(apply(cctfit$ppVDI,2,function(x) all(x==0)))){
cctfit$ppVDI[,i] <- NA
}
}
for(v in 1:dim(cctfit$ppVDI)[2]){
if(identical(unique(cctfit$ppVDI[,v]),c(NA,0)) || identical(unique(cctfit$ppVDI[,v]),c(0,NA))){cctfit$ppVDI[,v] <- NA}
}
}
rm(vdi, varvec)
rm(eigv,leigv,tmp,tmp2,tmp3,ind);
cctfit$checksrun <- TRUE
}
if(cctfit$whmodel=="LTRM"){
if(cctfit$V==1){
cctfit$V <- 1;
cctfit$BUGSoutput$sims.list$Om <- array(1,c(cctfit$BUGSoutput$n.sims,cctfit$n));
if(length(dim(cctfit$BUGSoutput$sims.list[["T"]])) < 3){
cctfit$BUGSoutput$sims.list$T <- array(cctfit$BUGSoutput$sims.list[["T"]][,],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))}
if(length(dim(cctfit$BUGSoutput$sims.list[["gam"]])) < 3){
cctfit$BUGSoutput$sims.list$gam <- array(cctfit$BUGSoutput$sims.list[["gam"]][,],c(cctfit$BUGSoutput$n.sims,cctfit$C-1,1))}
if(cctfit$itemdiff==TRUE){
if(length(dim(cctfit$BUGSoutput$sims.list$lam)) < 3){
cctfit$BUGSoutput$sims.list$lam <- array(cctfit$BUGSoutput$sims.list[["lam"]][,],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))}
}
}
if(cctfit$itemdiff==FALSE){
if(cctfit$V==1){
cctfit$BUGSoutput$sims.list$lam <- array(1,c(cctfit$BUGSoutput$n.sims,cctfit$m,1))
}
if(cctfit$V > 1){
cctfit$BUGSoutput$sims.list$lam <- array(1,c(cctfit$BUGSoutput$n.sims,cctfit$m,cctfit$V))
}
}
cctfit$ppY <- array(NA, c(cctfit$n,cctfit$m,subsetsize));
cctfit$ppX <- array(NA, c(cctfit$n,cctfit$m,subsetsize));
for(samp in indices){
tautmp <- (cctfit$BUGSoutput$sims.list[["E"]][samp,]*t(exp(cctfit$BUGSoutput$sims.list[["lam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])))^-2
cctfit$ppX[,,which(indices==samp)] <- matrix(
rnorm((cctfit$n*cctfit$m),t(cctfit$BUGSoutput$sims.list[["T"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]]),(tautmp^(-.5))),
cctfit$n,cctfit$m)
deltatmp <- cctfit$BUGSoutput$sims.list[["a"]][samp,]*t(cctfit$BUGSoutput$sims.list[["gam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])+cctfit$BUGSoutput$sims.list[["b"]][samp,]
rc <- which(cctfit$ppX[,,which(indices==samp)] < deltatmp[,1],arr.ind=TRUE)
cctfit$ppY[cbind(rc,array(which(indices==samp),dim(rc)[1]))] <- 1
for(c in 1:(cctfit$C-1)){
rc <- which(cctfit$ppX[,,which(indices==samp)] > deltatmp[,c],arr.ind=TRUE)
cctfit$ppY[cbind(rc,array(which(indices==samp),dim(rc)[1]))] <- (c+1) }
}
cctfit$mean$ppY <- rowMeans(cctfit$ppY[,,],dims=2)
if(cctfit$mval==TRUE){
for(i in 1:dim(cctfit$datob$thena)[1]){cctfit$data[cctfit$datob$thena[i,1],cctfit$datob$thena[i,2]] <- Mode(cctfit$ppY[cctfit$datob$thena[i,1],cctfit$datob$thena[i,2],])}
cctfit$MVest <- cbind(cctfit$datob$thena,cctfit$data[cctfit$datob$thena])
colnames(cctfit$MVest) <- c("Pers","Item","Resp")
}
leigv <- 12
eigv <- matrix(-1,dim(cctfit$ppY)[3],leigv)
options(warn=-3)
ind <- sample(indices,min(cctfit$BUGSoutput$n.sims,250,subsetsize)) #250 is the number of samples
tmp <- apply(cctfit$ppY,c(1,3),function(x) t(x))
if(polych==TRUE){
message(" note: utilizing polychoric correlations (time intensive)");
tmp2 <- apply(tmp[,,],3,function(x) suppressMessages(polychoric(x,global=FALSE))$rho);
tmp2 <- sapply(tmp2,function(x) array(x,c(cctfit$n,cctfit$n))) }else{
tmp2 <- apply(tmp[,,],3,function(x) cor(x))}
tmp3 <- array(tmp2,c(cctfit$n,cctfit$n,length(ind)))
for(i in 1:length(ind)){
suppressMessages(try(eigv[i,] <- fa(tmp3[,,i])$values[1:leigv],silent=TRUE))}
wch <- -which(eigv[,1]==-1);
if(length(wch)==0){cctfit$ppeig <- eigv}else{cctfit$ppeig <- eigv[-which(eigv[,1]==-1),]}
if(polych==TRUE){
cctfit$dateig <- suppressMessages(fa(polychoric(t(cctfit$data),global=FALSE)$rho)$values[1:leigv])}else{
cctfit$dateig <- suppressMessages(fa(cor(t(cctfit$data)))$values[1:leigv])
}
options(warn=0)
if(cctfit$V==1){
varvec <- matrix(NA,length(ind),dim(cctfit$ppY)[2])
vdi <- matrix(NA,dim(cctfit$ppY)[3],1)
for(i in 1:length(ind)){
varvec[i,] <- apply(cctfit$ppY[,,i],2,var)
}
vdi <- apply(varvec,1,var)
cctfit$ppVDI <- vdi
cctfit$datVDI <- var(apply(cctfit$data,2,var))
}
options(warn=-3)
if(cctfit$V>1){
varvec <- array(NA,c(length(ind),dim(cctfit$ppY)[2],cctfit$V))
vdi <- matrix(NA,dim(cctfit$ppY)[3],cctfit$V)
cctfit$datVDI <- array(-1,c(1,cctfit$V))
for(v in 1:cctfit$V){
for(i in 1:length(ind)){
if(sum(cctfit$BUGSoutput$sims.list$Om[ind[i],]==v)>1){
suppressMessages(try(varvec[i,,v] <- apply(cctfit$ppY[cctfit$BUGSoutput$sims.list$Om[ind[i],]==v,,i],2,var),silent=TRUE))
}else{
suppressMessages(try(varvec[i,,v] <- var(cctfit$ppY[cctfit$BUGSoutput$sims.list$Om[ind[i],]==v,,i]),silent=TRUE))
}
}
if(sum(cctfit$respmem==v)>1){
cctfit$datVDI[v] <- var(apply(cctfit$data[cctfit$respmem==v,],2,var))
}else{
cctfit$datVDI[v] <- var(cctfit$data[cctfit$respmem==v,])
}
}
vdi <- apply(varvec,c(1,3),var)
cctfit$ppVDI <- vdi
options(warn=0)
if(sum(apply(cctfit$ppVDI,2,function(x) all(x==0,na.rm=TRUE)))>=1){
for(i in which(apply(cctfit$ppVDI,2,function(x) all(x==0)))){
cctfit$ppVDI[,i] <- NA
}
}
for(v in 1:dim(cctfit$ppVDI)[2]){
if(identical(unique(cctfit$ppVDI[,v]),c(NA,0)) || identical(unique(cctfit$ppVDI[,v]),c(0,NA))){cctfit$ppVDI[,v] <- NA}
}
}
cctfit$polycor <- polych
rm(vdi, varvec)
rm(eigv,leigv,tmp,tmp2,tmp3,ind);
cctfit$checksrun <- TRUE
}
if(cctfit$whmodel=="CRM"){
invlogit <- function(x){x <- 1 / (1 + exp(-x)); return(x)}
logit <- function(x){x <- log(x/(1-x)); return(x)}
usetransform <- 0
version2 <- 1
if(cctfit$V==1){
cctfit$V <- 1;
cctfit$BUGSoutput$sims.list$Om <- array(1,c(cctfit$BUGSoutput$n.sims,cctfit$n));
if(length(dim(cctfit$BUGSoutput$sims.list[["T"]])) < 3){
cctfit$BUGSoutput$sims.list$T <- array(cctfit$BUGSoutput$sims.list[["T"]][,],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))}
if(cctfit$itemdiff==TRUE){
if(length(dim(cctfit$BUGSoutput$sims.list$lam)) < 3){
cctfit$BUGSoutput$sims.list$lam <- array(cctfit$BUGSoutput$sims.list[["lam"]][,],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))}
}
}
if(cctfit$itemdiff==FALSE){
if(cctfit$V==1){
cctfit$BUGSoutput$sims.list$lam <- array(0,c(cctfit$BUGSoutput$n.sims,cctfit$m,1))
}
if(cctfit$V > 1){
cctfit$BUGSoutput$sims.list$lam <- array(0,c(cctfit$BUGSoutput$n.sims,cctfit$m,cctfit$V))
}
}
cctfit$ppY <- array(NA, c(cctfit$n,cctfit$m,subsetsize));
cctfit$ppX <- array(NA, c(cctfit$n,cctfit$m,subsetsize));
for(samp in indices){
tautmp <- (cctfit$BUGSoutput$sims.list[["E"]][samp,]*t(exp(cctfit$BUGSoutput$sims.list[["lam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])))^-2
cctfit$ppY[,,which(indices==samp)] <- matrix(
rnorm((cctfit$n*cctfit$m),(cctfit$BUGSoutput$sims.list[["a"]][samp,]*t(cctfit$BUGSoutput$sims.list[["T"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]]))+matrix(rep(cctfit$BUGSoutput$sims.list[["b"]][samp,],cctfit$m),cctfit$n,cctfit$m),
(cctfit$BUGSoutput$sims.list[["a"]][samp,]*(tautmp^(-.5)))),
cctfit$n,cctfit$m)
cctfit$ppX[,,which(indices==samp)] <- matrix(
rnorm((cctfit$n*cctfit$m),t(cctfit$BUGSoutput$sims.list[["T"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]]),tautmp^(-.5)),c(cctfit$n,cctfit$m))
}
if(usetransform==1){
if(sum(cctfit$ppX < logit(.001)) > 0){cctfit$ppX[cctfit$ppX < logit(.001)] <- logit(.001)}
if(sum(cctfit$ppY < logit(.001)) > 0){cctfit$ppY[cctfit$ppY < logit(.001)] <- logit(.001)}
if(sum(cctfit$ppX > logit(.999)) > 0){cctfit$ppX[cctfit$ppX > logit(.999)] <- logit(.999)}
if(sum(cctfit$ppY > logit(.999)) > 0){cctfit$ppY[cctfit$ppY > logit(.999)] <- logit(.999)}
}
if(usetransform==2){
if(sum(cctfit$ppX < logit(.0001)) > 0){cctfit$ppX[cctfit$ppX < logit(.0001)] <- logit(.0001)}
if(sum(cctfit$ppY < logit(.0001)) > 0){cctfit$ppY[cctfit$ppY < logit(.0001)] <- logit(.0001)}
if(sum(cctfit$ppX > logit(.9999)) > 0){cctfit$ppX[cctfit$ppX > logit(.9999)] <- logit(.9999)}
if(sum(cctfit$ppY > logit(.9999)) > 0){cctfit$ppY[cctfit$ppY > logit(.9999)] <- logit(.9999)}
}
cctfit$mean$ppY <- rowMeans(cctfit$ppY[,,],dims=2)
if(cctfit$mval==TRUE){
for(i in 1:dim(cctfit$datob$thena)[1]){cctfit$data[cctfit$datob$thena[i,1],cctfit$datob$thena[i,2]] <- mean(cctfit$ppY[cctfit$datob$thena[i,1],cctfit$datob$thena[i,2],])}
cctfit$MVest <- cbind(cctfit$datob$thena,cctfit$data[cctfit$datob$thena])
colnames(cctfit$MVest) <- c("Pers","Item","Resp")
}
leigv <- 12
options(warn=-3)
ind <- indices #500 is the number of samples
eigv <- matrix(-1,length(ind),leigv)
tmp <- apply(cctfit$ppY,c(1,3),function(x) t(x))
tmp2 <- apply(tmp[,,],3,function(x) cor(x))
tmp3 <- array(tmp2,c(cctfit$n,cctfit$n,subsetsize))
for(i in 1:length(ind)){
suppressMessages(try(eigv[i,] <- fa(tmp3[,,i])$values[1:leigv],silent=TRUE))}
wch <- -which(eigv[,1]==-1);
if(length(wch)==0){cctfit$ppeig <- eigv}else{cctfit$ppeig <- eigv[-which(eigv[,1]==-1),]}
cctfit$dateig <- suppressMessages(fa(cor(t(cctfit$data)))$values[1:leigv])
options(warn=0)
#VDI Check
usetransform <- 0
if(usetransform==0){
invlogit <- function(x){return(x)}
}
if(cctfit$V==1){
vdi <- matrix(-1,dim(cctfit$ppY)[3],1)
varvec <- array(-1,c(length(ind),dim(cctfit$ppY)[2]))
for(i in 1:length(ind)){
varvec[i,] <- apply(invlogit(cctfit$ppY[,,i]),2,var)
if(version2==1){
varvec[i,] <- apply(invlogit(cctfit$ppX[,,i]),2,var)
}
}
vdi <- apply(varvec,1,var)
cctfit$ppVDI <- apply(varvec,1,var)
cctfit$datVDI <- var(apply(invlogit(cctfit$data),2,var))
if(version2==1){
cctfit$datVDI <- var(apply(invlogit((array(rep(1/cctfit$BUGSoutput$mean$a,times=cctfit$m),c(cctfit$n,cctfit$m))*cctfit$data)-matrix(rep(cctfit$BUGSoutput$mean$b,cctfit$m),cctfit$n,cctfit$m)),2,var))
}
}
if(cctfit$V>1){
varvec <- array(NA,c(length(ind),dim(cctfit$ppY)[2],cctfit$V))
cctfit$datVDI <- array(NA,c(1,cctfit$V))
options(warn=-3)
datest <- t(sapply(1:cctfit$n,function(i) ((1/mean(cctfit$BUGSoutput$sims.list$a[cctfit$BUGSoutput$sims.list$Om[,i]==cctfit$respmem[i],i]))*
cctfit$data[i,])-(mean(cctfit$BUGSoutput$sims.list$b[cctfit$BUGSoutput$sims.list$Om[,i]==cctfit$respmem[i],i]))
))
for(v in 1:cctfit$V){
for(i in 1:length(ind)){
if(sum(cctfit$BUGSoutput$sims.list$Om[ind[i],]==v)>1){
suppressMessages(try(varvec[i,,v] <- apply(invlogit(cctfit$ppY[cctfit$BUGSoutput$sims.list$Om[ind[i],]==v,,i]),2,var),silent=TRUE))
if(version2==1){
suppressMessages(try(varvec[i,,v] <- apply(invlogit(cctfit$ppX[cctfit$BUGSoutput$sims.list$Om[ind[i],]==v,,i]),2,var),silent=TRUE))
}
}else{
suppressMessages(try(varvec[i,,v] <- var(invlogit(cctfit$ppY[cctfit$BUGSoutput$sims.list$Om[ind[i],]==v,,i])),silent=TRUE))
if(version2==1){
suppressMessages(try(varvec[i,,v] <- var(invlogit(cctfit$ppX[cctfit$BUGSoutput$sims.list$Om[ind[i],]==v,,i])),silent=TRUE))
}
}
}
if(sum(cctfit$respmem==v)>1){
cctfit$datVDI[v] <- var(apply(datest[cctfit$respmem==v ,],2,var))
}else{
cctfit$datVDI[v] <- var(datest[cctfit$respmem==v ,])
}
}
options(warn=0)
vdi <- apply(varvec,c(1,3),function(x) var(x,na.rm=TRUE))
cctfit$ppVDI <- vdi
if(sum(apply(cctfit$ppVDI,2,function(x) all(x==0,na.rm=TRUE)))>=1){
for(i in which(apply(cctfit$ppVDI,2,function(x) all(x==0)))){
cctfit$ppVDI[,i] <- NA
}
}
for(v in 1:dim(cctfit$ppVDI)[2]){
if(identical(unique(cctfit$ppVDI[,v]),c(NA,0)) || identical(unique(cctfit$ppVDI[,v]),c(0,NA))){cctfit$ppVDI[,v] <- NA}
}
}
if(usetransform==0){
invlogit <- function(x){x <- 1 / (1 + exp(-x)); return(x)}
}
rm(vdi, varvec)
rm(eigv,leigv,tmp,tmp2,tmp3,ind);
cctfit$checksrun <- TRUE
}
}else{
#########################################
##### Calculate Posterior Predictive Data From All Samples (memory intensive)
#########################################
if(cctfit$whmodel=="GCM"){
if(cctfit$V==1){
cctfit$V <- 1;
cctfit$BUGSoutput$sims.list$Om <- array(1,c(cctfit$BUGSoutput$n.sims,cctfit$n));
if(length(dim(cctfit$BUGSoutput$sims.list$Z)) < 3){
cctfit$BUGSoutput$sims.list$Z <- array(cctfit$BUGSoutput$sims.list[["Z"]][,],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))}
if(cctfit$itemdiff==TRUE){
if(length(dim(cctfit$BUGSoutput$sims.list$lam)) < 3){
cctfit$BUGSoutput$sims.list$lam <- array(cctfit$BUGSoutput$sims.list[["lam"]][,],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))}
}
}
if(cctfit$itemdiff==FALSE){
if(cctfit$V==1){
cctfit$BUGSoutput$sims.list$lam <- array(.5,c(cctfit$BUGSoutput$n.sims,cctfit$m,1))
}
if(cctfit$V > 1){
cctfit$BUGSoutput$sims.list$lam <- array(.5,c(cctfit$BUGSoutput$n.sims,cctfit$m,cctfit$V))
} }
cctfit$ppY <- array(NA, c(cctfit$n,cctfit$m,cctfit$BUGSoutput$n.sims));
cctfit$BUGSoutput$sims.list$tau <- array(-1,c(cctfit$BUGSoutput$n.sims,cctfit$n,cctfit$m))
for(samp in 1:cctfit$BUGSoutput$n.sims){
cctfit$BUGSoutput$sims.list[["tau"]][samp,,] <- (cctfit$BUGSoutput$sims.list[["th"]][samp,]*t(1-cctfit$BUGSoutput$sims.list[["lam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])) /
( (cctfit$BUGSoutput$sims.list[["th"]][samp,]*t(1-cctfit$BUGSoutput$sims.list[["lam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])) +
((1-cctfit$BUGSoutput$sims.list[["th"]][samp,])*t(cctfit$BUGSoutput$sims.list[["lam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])) )
cctfit$ppY[,,samp] <- matrix(rbinom((cctfit$n*cctfit$m),1,
(t(cctfit$BUGSoutput$sims.list[["Z"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])*cctfit$BUGSoutput$sims.list[["tau"]][samp,,] ) +
t(t(1-cctfit$BUGSoutput$sims.list[["tau"]][samp,,])%*%diag(cctfit$BUGSoutput$sims.list[["g"]][samp,])) ), cctfit$n,cctfit$m)
}
if(cctfit$mval==TRUE){
for(i in 1:dim(cctfit$datob$thena)[1]){cctfit$data[cctfit$datob$thena[i,1],cctfit$datob$thena[i,2]] <- Mode(cctfit$ppY[cctfit$datob$thena[i,1],cctfit$datob$thena[i,2],])}
cctfit$MVest <- cbind(cctfit$datob$thena,cctfit$data[cctfit$datob$thena])
colnames(cctfit$MVest) <- c("Pers","Item","Resp")
}
cctfit$mean$ppY <- rowMeans(cctfit$ppY[,,],dims=2)
leigv <- 12
options(warn=-3)
ind <- sample(cctfit$BUGSoutput$n.sims,min(500,cctfit$BUGSoutput$n.sims)) #500 is the number of samples
eigv <- matrix(-1,length(ind),leigv)
tmp <- apply(cctfit$ppY[,,ind],c(1,3),function(x) t(x))
tmp2 <- apply(tmp[,,],3,function(x) cor(x))
tmp3 <- array(tmp2,c(cctfit$n,cctfit$n,cctfit$BUGSoutput$n.sims))
for(i in 1:length(ind)){
suppressMessages(try(eigv[i,] <- fa(tmp3[,,i])$values[1:leigv],silent=TRUE))}
wch <- -which(eigv[,1]==-1);
if(length(wch)==0){cctfit$ppeig <- eigv}else{cctfit$ppeig <- eigv[-which(eigv[,1]==-1),]}
cctfit$dateig <- suppressMessages(fa(cor(t(cctfit$data)))$values[1:leigv])
options(warn=0)
if(cctfit$V==1){
varvec <- matrix(-1,length(ind),dim(cctfit$ppY)[2])
vdi <- matrix(-1,dim(cctfit$ppY)[3],1)
for(i in 1:length(ind)){
varvec[i,] <- apply(cctfit$ppY[,,i],2,var)
}
vdi <- apply(varvec,1,var)
cctfit$ppVDI <- vdi
cctfit$datVDI <- var(apply(cctfit$data,2,var))
}
options(warn=-3)
if(cctfit$V>1){
varvec <- array(NA,c(length(ind),dim(cctfit$ppY)[2],cctfit$V))
vdi <- matrix(NA,dim(cctfit$ppY)[3],cctfit$V)
cctfit$datVDI <- array(-1,c(1,cctfit$V))
for(v in 1:cctfit$V){
for(i in 1:length(ind)){
if(sum(cctfit$BUGSoutput$sims.list$Om[ind[i],]==v)>1){
suppressMessages(try(varvec[i,,v] <- apply(cctfit$ppY[cctfit$BUGSoutput$sims.list$Om[ind[i],]==v,,i],2,var),silent=TRUE))
}else{
suppressMessages(try(varvec[i,,v] <- var(cctfit$ppY[cctfit$BUGSoutput$sims.list$Om[ind[i],]==v,,i]),silent=TRUE))
}
}
if(sum(cctfit$respmem==v)>1){
cctfit$datVDI[v] <- var(apply(cctfit$data[cctfit$respmem==v,],2,var))
}else{
cctfit$datVDI[v] <- var(cctfit$data[cctfit$respmem==v,])
}
}
vdi <- apply(varvec,c(1,3),var)
cctfit$ppVDI <- vdi
options(warn=0)
if(sum(apply(cctfit$ppVDI,2,function(x) all(x==0,na.rm=TRUE)))>=1){
for(i in which(apply(cctfit$ppVDI,2,function(x) all(x==0)))){
cctfit$ppVDI[,i] <- NA
}
}
for(v in 1:dim(cctfit$ppVDI)[2]){
if(identical(unique(cctfit$ppVDI[,v]),c(NA,0)) || identical(unique(cctfit$ppVDI[,v]),c(0,NA))){cctfit$ppVDI[,v] <- NA}
}
}
rm(vdi, varvec)
rm(eigv,leigv,tmp,tmp2,tmp3,ind);
cctfit$checksrun <- TRUE
}
if(cctfit$whmodel=="LTRM"){
if(cctfit$V==1){
cctfit$V <- 1;
cctfit$BUGSoutput$sims.list$Om <- array(1,c(cctfit$BUGSoutput$n.sims,cctfit$n));
if(length(dim(cctfit$BUGSoutput$sims.list[["T"]])) < 3){
cctfit$BUGSoutput$sims.list$T <- array(cctfit$BUGSoutput$sims.list[["T"]][,],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))}
if(length(dim(cctfit$BUGSoutput$sims.list[["gam"]])) < 3){
cctfit$BUGSoutput$sims.list$gam <- array(cctfit$BUGSoutput$sims.list[["gam"]][,],c(cctfit$BUGSoutput$n.sims,cctfit$C-1,1))}
if(cctfit$itemdiff==TRUE){
if(length(dim(cctfit$BUGSoutput$sims.list$lam)) < 3){
cctfit$BUGSoutput$sims.list$lam <- array(cctfit$BUGSoutput$sims.list[["lam"]][,],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))}
}
}
if(cctfit$itemdiff==FALSE){
if(cctfit$V==1){
cctfit$BUGSoutput$sims.list$lam <- array(1,c(cctfit$BUGSoutput$n.sims,cctfit$m,1))
}
if(cctfit$V > 1){
cctfit$BUGSoutput$sims.list$lam <- array(1,c(cctfit$BUGSoutput$n.sims,cctfit$m,cctfit$V))
}
}
cctfit$ppY <- array(NA, c(cctfit$n,cctfit$m,cctfit$BUGSoutput$n.sims));
cctfit$ppX <- array(NA, c(cctfit$n,cctfit$m,cctfit$BUGSoutput$n.sims));
cctfit$ppdelta <- array(NA, c(cctfit$n,cctfit$C-1,cctfit$BUGSoutput$n.sims))
cctfit$BUGSoutput$sims.list$tau <- array(-1,c(cctfit$BUGSoutput$n.sims,cctfit$n,cctfit$m))
for(samp in 1:cctfit$BUGSoutput$n.sims){
cctfit$BUGSoutput$sims.list[["tau"]][samp,,] <- cctfit$BUGSoutput$sims.list[["E"]][samp,]*t(1/cctfit$BUGSoutput$sims.list[["lam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])
cctfit$ppX[,,samp] <- matrix(
rnorm((cctfit$n*cctfit$m),t(cctfit$BUGSoutput$sims.list[["T"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]]),(cctfit$BUGSoutput$sims.list[["tau"]][samp,,]^(-.5))),
cctfit$n,cctfit$m)
cctfit$ppdelta[,,samp] <- cctfit$BUGSoutput$sims.list[["a"]][samp,]*t(cctfit$BUGSoutput$sims.list[["gam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])+cctfit$BUGSoutput$sims.list[["b"]][samp,]
rc <- which(cctfit$ppX[,,samp] < cctfit$ppdelta[,1,samp],arr.ind=TRUE)
cctfit$ppY[cbind(rc,array(samp,dim(rc)[1]))] <- 1
for(c in 1:(cctfit$C-1)){
rc <- which(cctfit$ppX[,,samp] > cctfit$ppdelta[,c,samp],arr.ind=TRUE)
cctfit$ppY[cbind(rc,array(samp,dim(rc)[1]))] <- (c+1) }
}
cctfit$mean$ppY <- rowMeans(cctfit$ppY[,,],dims=2)
if(cctfit$mval==TRUE){
for(i in 1:dim(cctfit$datob$thena)[1]){cctfit$data[cctfit$datob$thena[i,1],cctfit$datob$thena[i,2]] <- Mode(cctfit$ppY[cctfit$datob$thena[i,1],cctfit$datob$thena[i,2],])}
cctfit$MVest <- cbind(cctfit$datob$thena,cctfit$data[cctfit$datob$thena])
colnames(cctfit$MVest) <- c("Pers","Item","Resp")
}
leigv <- 12
eigv <- matrix(-1,dim(cctfit$ppY)[3],leigv)
options(warn=-3)
ind <- sample(cctfit$BUGSoutput$n.sims,min(250,cctfit$BUGSoutput$n.sims)) #250 is the number of samples
tmp <- apply(cctfit$ppY[,,ind],c(1,3),function(x) t(x))
if(polych==TRUE){tmp2 <- apply(tmp[,,],3,function(x) suppressMessages(polychoric(x,global=FALSE))$rho)}else{
tmp2 <- apply(tmp[,,],3,function(x) cor(x))
}
tmp3 <- array(tmp2,c(cctfit$n,cctfit$n,length(ind)))
for(i in 1:length(ind)){
suppressMessages(try(eigv[i,] <- fa(tmp3[,,i])$values[1:leigv],silent=TRUE))}
wch <- -which(eigv[,1]==-1);
if(length(wch)==0){cctfit$ppeig <- eigv}else{cctfit$ppeig <- eigv[-which(eigv[,1]==-1),]}
if(polych==TRUE){
cctfit$dateig <- suppressMessages(fa(polychoric(t(cctfit$data),global=FALSE)$rho)$values[1:leigv])}else{
cctfit$dateig <- suppressMessages(fa(cor(t(cctfit$data)))$values[1:leigv])
}
options(warn=0)
if(cctfit$V==1){
varvec <- matrix(NA,length(ind),dim(cctfit$ppY)[2])
vdi <- matrix(NA,dim(cctfit$ppY)[3],1)
for(i in 1:length(ind)){
varvec[i,] <- apply(cctfit$ppY[,,i],2,var)
}
vdi <- apply(varvec,1,var)
cctfit$ppVDI <- vdi
cctfit$datVDI <- var(apply(cctfit$data,2,var))
}
options(warn=-3)
if(cctfit$V>1){
varvec <- array(NA,c(length(ind),dim(cctfit$ppY)[2],cctfit$V))
vdi <- matrix(NA,dim(cctfit$ppY)[3],cctfit$V)
cctfit$datVDI <- array(-1,c(1,cctfit$V))
for(v in 1:cctfit$V){
for(i in 1:length(ind)){
if(sum(cctfit$BUGSoutput$sims.list$Om[ind[i],]==v)>1){
suppressMessages(try(varvec[i,,v] <- apply(cctfit$ppY[cctfit$BUGSoutput$sims.list$Om[ind[i],]==v,,i],2,var),silent=TRUE))
}else{
suppressMessages(try(varvec[i,,v] <- var(cctfit$ppY[cctfit$BUGSoutput$sims.list$Om[ind[i],]==v,,i]),silent=TRUE))
}
}
if(sum(cctfit$respmem==v)>1){
cctfit$datVDI[v] <- var(apply(cctfit$data[cctfit$respmem==v,],2,var))
}else{
cctfit$datVDI[v] <- var(cctfit$data[cctfit$respmem==v,])
}
}
vdi <- apply(varvec,c(1,3),var)
cctfit$ppVDI <- vdi
options(warn=0)
if(sum(apply(cctfit$ppVDI,2,function(x) all(x==0,na.rm=TRUE)))>=1){
for(i in which(apply(cctfit$ppVDI,2,function(x) all(x==0)))){
cctfit$ppVDI[,i] <- NA
}
}
for(v in 1:dim(cctfit$ppVDI)[2]){
if(identical(unique(cctfit$ppVDI[,v]),c(NA,0)) || identical(unique(cctfit$ppVDI[,v]),c(0,NA))){cctfit$ppVDI[,v] <- NA}
}
}
cctfit$polycor <- polych
rm(vdi, varvec)
rm(eigv,leigv,tmp,tmp2,tmp3,ind);
cctfit$checksrun <- TRUE
}
if(cctfit$whmodel=="CRM"){
invlogit <- function(x){x <- 1 / (1 + exp(-x)); return(x)}
logit <- function(x){x <- log(x/(1-x)); return(x)}
if(cctfit$V==1){
cctfit$V <- 1;
cctfit$BUGSoutput$sims.list$Om <- array(1,c(cctfit$BUGSoutput$n.sims,cctfit$n));
if(length(dim(cctfit$BUGSoutput$sims.list[["T"]])) < 3){
cctfit$BUGSoutput$sims.list$T <- array(cctfit$BUGSoutput$sims.list[["T"]][,],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))}
if(cctfit$itemdiff==TRUE){
if(length(dim(cctfit$BUGSoutput$sims.list$lam)) < 3){
cctfit$BUGSoutput$sims.list$lam <- array(cctfit$BUGSoutput$sims.list[["lam"]][,],c(cctfit$BUGSoutput$n.sims,cctfit$m,1))}
}
}
if(cctfit$itemdiff==FALSE){
if(cctfit$V==1){
cctfit$BUGSoutput$sims.list$lam <- array(1,c(cctfit$BUGSoutput$n.sims,cctfit$m,1))
}
if(cctfit$V > 1){
cctfit$BUGSoutput$sims.list$lam <- array(1,c(cctfit$BUGSoutput$n.sims,cctfit$m,cctfit$V))
}
}
cctfit$ppY <- array(NA, c(cctfit$n,cctfit$m,cctfit$BUGSoutput$n.sims));
cctfit$ppX <- array(NA, c(cctfit$n,cctfit$m,cctfit$BUGSoutput$n.sims));
cctfit$BUGSoutput$sims.list$tau <- array(-1,c(cctfit$BUGSoutput$n.sims,cctfit$n,cctfit$m))
for(samp in 1:cctfit$BUGSoutput$n.sims){
cctfit$BUGSoutput$sims.list[["tau"]][samp,,] <- cctfit$BUGSoutput$sims.list[["E"]][samp,]*t(1/cctfit$BUGSoutput$sims.list[["lam"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]])
cctfit$ppY[,,samp] <- matrix(
rnorm((cctfit$n*cctfit$m),(cctfit$BUGSoutput$sims.list[["a"]][samp,]*t(cctfit$BUGSoutput$sims.list[["T"]][samp,,cctfit$BUGSoutput$sims.list[["Om"]][samp,]]))+matrix(rep(cctfit$BUGSoutput$sims.list[["b"]][samp,],cctfit$m),cctfit$n,cctfit$m),
(cctfit$BUGSoutput$sims.list[["a"]][samp,]*(cctfit$BUGSoutput$sims.list[["tau"]][samp,,]^(-.5)))),
cctfit$n,cctfit$m)
cctfit$ppX[,,samp] <- (array(rep(1/cctfit$BUGSoutput$sims.list[["a"]][samp,],times=cctfit$m),c(cctfit$n,cctfit$m))*(cctfit$ppY[,,samp]))-matrix(rep(cctfit$BUGSoutput$sims.list[["b"]][samp,],times=cctfit$m),c(cctfit$n,cctfit$m))
}
cctfit$mean$ppY <- rowMeans(cctfit$ppY[,,],dims=2)
if(cctfit$mval==TRUE){
for(i in 1:dim(cctfit$datob$thena)[1]){cctfit$data[cctfit$datob$thena[i,1],cctfit$datob$thena[i,2]] <- mean(cctfit$ppY[cctfit$datob$thena[i,1],cctfit$datob$thena[i,2],])}
cctfit$MVest <- cbind(cctfit$datob$thena,cctfit$data[cctfit$datob$thena])
colnames(cctfit$MVest) <- c("Pers","Item","Resp")
}
leigv <- 12
options(warn=-3)
ind <- sample(cctfit$BUGSoutput$n.sims,min(500,cctfit$BUGSoutput$n.sims)) #500 is the number of samples
eigv <- matrix(-1,length(ind),leigv)
tmp <- apply(cctfit$ppY[,,ind],c(1,3),function(x) t(x))
tmp2 <- apply(tmp[,,],3,function(x) cor(x))
tmp3 <- array(tmp2,c(cctfit$n,cctfit$n,cctfit$BUGSoutput$n.sims))
for(i in 1:length(ind)){
suppressMessages(try(eigv[i,] <- fa(tmp3[,,i])$values[1:leigv],silent=TRUE))}
wch <- -which(eigv[,1]==-1);
if(length(wch)==0){cctfit$ppeig <- eigv}else{cctfit$ppeig <- eigv[-which(eigv[,1]==-1),]}
cctfit$dateig <- suppressMessages(fa(cor(t(cctfit$data)))$values[1:leigv])
options(warn=0)
if(cctfit$V==1){
vdi <- matrix(-1,dim(cctfit$ppY)[3],1)
varvec <- array(-1,c(length(ind),dim(cctfit$ppY)[2]))
for(i in 1:length(ind)){
#varvec[i,] <- apply(invlogit(cctfit$ppY[,,i]),2,var)
varvec[i,] <- apply(invlogit(cctfit$ppY[,,i]),2,var)
}
vdi <- apply(varvec,1,var)
cctfit$ppVDI <- apply(varvec,1,var)
#cctfit$datVDI <- var(apply(invlogit(cctfit$data),2,var))
cctfit$datVDI <- var(apply(invlogit(cctfit$data),2,var))
}
if(cctfit$V>1){
varvec <- array(NA,c(length(ind),dim(cctfit$ppY)[2],cctfit$V))
cctfit$datVDI <- array(NA,c(1,cctfit$V))
options(warn=-3)
for(v in 1:cctfit$V){
for(i in 1:length(ind)){
if(sum(cctfit$BUGSoutput$sims.list$Om[ind[i],]==v)>1){
#suppressMessages(try(varvec[i,,v] <- apply(invlogit(cctfit$ppY[cctfit$BUGSoutput$sims.list$Om[ind[i],]==v,,i]),2,var),silent=TRUE))
suppressMessages(try(varvec[i,,v] <- apply(invlogit(cctfit$ppY[cctfit$BUGSoutput$sims.list$Om[ind[i],]==v,,i]),2,var),silent=TRUE))
}else{
#suppressMessages(try(varvec[i,,v] <- var(invlogit(cctfit$ppY[cctfit$BUGSoutput$sims.list$Om[ind[i],]==v,,i])),silent=TRUE))
suppressMessages(try(varvec[i,,v] <- var(invlogit(cctfit$ppY[cctfit$BUGSoutput$sims.list$Om[ind[i],]==v,,i])),silent=TRUE))
}
}
if(sum(cctfit$respmem==v)>1){
#cctfit$datVDI[v] <- var(apply(invlogit(cctfit$data[cctfit$respmem==v,]),2,var))
cctfit$datVDI[v] <- var(apply(invlogit(cctfit$data[cctfit$respmem==v,]),2,var))
}else{
#cctfit$datVDI[v] <- var(invlogit(cctfit$data[cctfit$respmem==v,]))
cctfit$datVDI[v] <- var(invlogit(cctfit$data[cctfit$respmem==v,]))
}
}
options(warn=0)
vdi <- apply(varvec,c(1,3),function(x) var(x,na.rm=TRUE))
cctfit$ppVDI <- vdi
if(sum(apply(cctfit$ppVDI,2,function(x) all(x==0,na.rm=TRUE)))>=1){
for(i in which(apply(cctfit$ppVDI,2,function(x) all(x==0)))){
cctfit$ppVDI[,i] <- NA
}
}
for(v in 1:dim(cctfit$ppVDI)[2]){
if(identical(unique(cctfit$ppVDI[,v]),c(NA,0)) || identical(unique(cctfit$ppVDI[,v]),c(0,NA))){cctfit$ppVDI[,v] <- NA}
}
}
rm(vdi, varvec)
rm(eigv,leigv,tmp,tmp2,tmp3,ind);
cctfit$checksrun <- TRUE
}
}
if(cctfit$mval==TRUE && gui==FALSE){
message(paste("\n ...Use function 'cctmvest()' to view model estimates for missing data",sep=""))
# message(paste("\n ...Use '",guidat$varname,"$MVest' to view posterior predictive estimates for missing data\n ... '",guidat$varname,"$data' to view the full data matrix \n ... '",guidat$varname,"$datamiss' to view the matrix with missing values",sep=""))
}
message("\n ...Posterior predictive checks complete")
}
if(doplot == TRUE){
if(saveplots==1){jpeg(file.path(gsub(".Rdata","ppc.jpg",savedir)),width = 6, height = 3, units = "in", pointsize = 12,quality=100,res=400)}
if(saveplots==2){postscript(file=file.path(gsub(".Rdata","ppc.eps",savedir)), onefile=FALSE, horizontal=FALSE, width = 6, height = 6, paper="special", family="Times")}
par(oma=c(0,0,0,0),mar=c(4,4,3,1),mgp=c(2.25,.75,0),mfrow=c(1,2))
plot(NA, main="Culture Number Check",xlim=c(1,min(cctfit$n,8)),ylim=c(min(cctfit$dateig[min(cctfit$n,8)],cctfit$ppeig[,min(cctfit$n,8)]),ceiling(max(cctfit$ppeig[,1],cctfit$dateig[1]))),xlab="Eigenvalue",ylab="Value",las=1,pch=21,type="l",col="black")
for(i in 1:dim(cctfit$ppeig)[1]){
points(cctfit$ppeig[i,],col="grey",type="l") }
points(cctfit$dateig,col="black",type="l")
if(cctfit$V>1){
cctfit$VDIperc <- array(NA,c(cctfit$V,1))
color <- "black"; color2 <- "black"
color <- c("black","dark grey","light grey",rainbow(cctfit$V))
minmax <- array(NA,c(cctfit$V,2))
minmax2 <- array(NA,c(cctfit$V,1))
for(v in which(sort(which(!apply(cctfit$ppVDI,2,function(x) all(is.na(x))))) %in% sort(unique(cctfit$respmem)))){
tmpdist <-ecdf(cctfit$ppVDI[,v])
minmax[v,] <- c(quantile(tmpdist, .005),quantile(tmpdist, .995))
minmax2[v] <- max(density(cctfit$ppVDI[,v],na.rm=TRUE)$y)
}
plot(NA,main="Item Difficulty Check", xlab= "VDI", xlim=c(min(cctfit$datVDI[!is.na(cctfit$datVDI)],minmax,na.rm=TRUE),max(cctfit$datVDI[!is.na(cctfit$datVDI)],minmax,na.rm=TRUE)), ylim=c(0,max(minmax2,na.rm=TRUE)), ylab="Value",las=1,col=color[1],type="l")
for(v in which(sort(which(!apply(cctfit$ppVDI,2,function(x) all(is.na(x))))) %in% sort(unique(cctfit$respmem)))){
tmpdist <-ecdf(cctfit$ppVDI[,v])
points(density(cctfit$ppVDI[,v],na.rm=TRUE),lwd=1.5,main="Item Difficulty Check", xlab= "VDI", xlim=c(min(cctfit$datVDI[!is.na(cctfit$datVDI)],quantile(tmpdist, .005),na.rm=TRUE),max(cctfit$datVDI[!is.na(cctfit$datVDI)],quantile(tmpdist, .995),na.rm=TRUE)), ylim=c(0,max(density(cctfit$ppVDI,na.rm=TRUE)$y,na.rm=TRUE)), ylab="Value",las=1,col=color[v],type="l")
if(saveplots!=1 && saveplots!=2){print(paste("VDI Culture ",v," : ",round(100*tmpdist(cctfit$datVDI[,v]),digits=2)," percentile"),col=color2)
cctfit$VDIperc[v,1] <- round(100*tmpdist(cctfit$datVDI[,v]),digits=2)
}
}
rm(tmpdist)
segments(cctfit$datVDI,0,cctfit$datVDI,par("usr")[4],col=color,lwd=1.5)
}else{
color <- "black"; color2 <- "black"
tmpdist <-ecdf(cctfit$ppVDI)
plot(density(cctfit$ppVDI,na.rm=TRUE),main="Item Difficulty Check", xlab= "VDI", xlim=c(min(cctfit$datVDI,quantile(tmpdist, .005)),max(cctfit$datVDI,quantile(tmpdist, .995))), ylim=c(0,max(density(cctfit$ppVDI,na.rm=TRUE)$y,na.rm=TRUE)), ylab="Value",las=1,col=color2)
segments(cctfit$datVDI,0,cctfit$datVDI,par("usr")[4],col=color2)
text(.7*par("usr")[2],.7*par("usr")[4],labels=paste(round(100*tmpdist(cctfit$datVDI),digits=2)," percentile"),col=color2)
if(saveplots!=1 && saveplots!=2){cctfit$VDIperc <- round(100*tmpdist(cctfit$datVDI),digits=2)
print(paste("VDI Culture 1 : ",round(100*tmpdist(cctfit$datVDI),digits=2)," percentile"),col=color2)
}; rm(tmpdist)
}
rm(color,color2)
if(saveplots==1 || saveplots==2){dev.off()}
saveplots <- 0
}
return(cctfit)
}
#######################
#Accessor Function for Missing Value Estimates
#######################
cctmvest <- function(cctfit){
if(cctfit$mval==FALSE){
message("\n ...Data has no missing values to estimate",sep="")
return()
}
if(cctfit$checksrun==FALSE){
message("\n ...Please use 'cctppc()' first",sep="")
}else{
return(cctfit$MVest)
}
}
#######################
#Accessor Function for the Factor Coefficients
#######################
cctfac <- function(dat){
return(print(dat$factors,digits=2))
}
#######################
#Accessor Function for the Subject Parameters
#######################
cctsubj <- function(cctfit){
return(print(cctfit$subj,digits=2))
}
#######################
#Accessor Function for the Credible Intervals of the Subject Parameters
#######################
cctsubjhdi <- function(cctfit){
return(print(cctfit$subjhdi,digits=2))
}
#######################
#Accessor Function for the Item Parameters
#######################
cctitem <- function(cctfit){
return(print(cctfit$item,digits=2))
}
#######################
#Accessor Function for the Credible Intervals of the Item Parameters
#######################
cctitemhdi <- function(cctfit){
return(print(cctfit$itemhdi,digits=2))
}
#######################
#Accessor Function for the Category Boundary Parameters and Credible Intervals (LTRM only)
#######################
cctcat <- function(cctfit){
return(print(cctfit$cat,digits=2))
}
#######################
#Accessor Function for Cluster Memberships
#######################
cctmemb <- function(cctfit){
return(cctfit$respmem)
}
#######################
#cctfit Summary Function
#######################
cctsum <- function(cctfit){
tmp <- c(" ","X")
if(cctfit$checksrun==FALSE){
cctfit$sumtab <- t(array(
c(paste("Data:"),cctfit$datob$datatype,paste("Model:"),cctfit$whmodel,
paste("Respondents:"),cctfit$n,paste("Items:"),cctfit$m,
paste("Observations: "),paste(100*round(cctfit$datob$nobs / (cctfit$datob$nobs+cctfit$datob$nmiss),2),"%",sep=""),paste("Missing Values: "),paste(100*round(cctfit$datob$nmiss / (cctfit$datob$nobs+cctfit$datob$nmiss),2),"%",sep=""),
paste("Cultures:"),cctfit$V,paste("Item Difficulty"),paste("[",tmp[cctfit$itemdiff+1],"]",sep=""),
paste("DIC:"),round(cctfit$BUGSoutput$DIC,2),paste("pD:"),round(cctfit$BUGSoutput$pD,2),
paste("Rhats > 1.10:"),paste(cctfit$Rhat$above110,"/",cctfit$Rhat$ncp),paste("Rhats > 1.05:"),paste(cctfit$Rhat$above105,"/",cctfit$Rhat$ncp),
paste("Samples:"),cctfit$BUGSoutput$n.iter,paste("Chains:"),cctfit$BUGSoutput$n.chains,
paste("Burn-in:"),cctfit$BUGSoutput$n.burnin,paste("Thinning:"),cctfit$BUGSoutput$n.thin
#,paste("PPC Calculated"),paste("[",tmp[cctfit$checksrun+1],"]",sep=""),paste(""),paste("")
),c(4,8)))
}else{
cctfit$sumtab <- t(array(
c(paste("Data:"),cctfit$datob$datatype,paste("Model:"),cctfit$whmodel,
paste("Respondents:"),cctfit$n,paste("Items:"),cctfit$m,
paste("Observations: "),paste(100*round(cctfit$datob$nobs / (cctfit$datob$nobs+cctfit$datob$nmiss),2),"%",sep=""),paste("Missing Values: "),paste(100*round(cctfit$datob$nmiss / (cctfit$datob$nobs+cctfit$datob$nmiss),2),"%",sep=""),
paste("Cultures:"),cctfit$V,paste("Item Difficulty"),paste("[",tmp[cctfit$itemdiff+1],"]",sep=""),
paste("DIC:"),round(cctfit$BUGSoutput$DIC,2),paste("pD:"),round(cctfit$BUGSoutput$pD,2),
paste("Rhats > 1.10:"),paste(cctfit$Rhat$above110,"/",cctfit$Rhat$ncp),paste("Rhats > 1.05:"),paste(cctfit$Rhat$above105,"/",cctfit$Rhat$ncp),
paste("PPC Calculated"),paste("[",tmp[cctfit$checksrun+1],"]",sep=""),paste("VDI"),paste(round(c(cctfit$VDI)),collapse=', '),
paste("Samples:"),cctfit$BUGSoutput$n.iter,paste("Chains:"),cctfit$BUGSoutput$n.chains,
paste("Burn-in:"),cctfit$BUGSoutput$n.burnin,paste("Thinning:"),cctfit$BUGSoutput$n.thin
),c(4,9)))
}
rownames(cctfit$sumtab) <- rep("", nrow(cctfit$sumtab))
colnames(cctfit$sumtab) <- rep("", ncol(cctfit$sumtab))
print(cctfit$sumtab,quote=F)
}
######################
#Function for the 'Export Results' Button
#- Prompts user where to save, and the filename to use
#- Saves the inference results (cctfit) to an .Rdata file of that filename
#- Saves .eps and .jpeg plots of the scree plot, plot results button, and run checks button
# if the checks were calculated (the checks button was pressed)
######################
exportfuncbutton <- function() {
guidat <- get("guidat", pkg_globals)
exportfunc(cctfit=get(guidat$varname, pkg_globals),gui=TRUE,guidat=guidat)
}
exportfunc <- function(cctfit,filename=paste(getwd(),"/CCTpackdata.Rdata",sep=""),gui=FALSE,guidat=NULL) {
if(gui==TRUE){
savedir <- tclvalue(tkgetSaveFile(initialfile = "CCTpackdata.Rdata",filetypes = "{{Rdata Files} {.Rdata}}"))
if (!nchar(savedir)) {
message("\n ...Export cancelled\n")
return()}
}else{
savedir <- filename
}
if (savedir=="" && !file.exists("CCTpack")){dir.create(file.path(getwd(), "CCTpack"));
savedir <- file.path(getwd(),"CCTpack","CCTpackdata.Rdata")
}
if(substr(savedir, nchar(savedir)-6+1, nchar(savedir))!=".Rdata"){
savedir <- paste(savedir,".Rdata",sep="")
}
message("\n ...Exporting results")
write.csv(cctfit$BUGSoutput$summary,file.path(gsub(".Rdata","posterior.csv",savedir)))
save(cctfit,file=file.path(savedir))
if(cctfit$checksrun==TRUE){
ppcfunc(cctfit=cctfit,saveplots=1,savedir); ppcfunc(cctfit=cctfit,saveplots=2,savedir)}
plotresultsfunc(cctfit=cctfit,saveplots=1,savedir); plotresultsfunc(cctfit=cctfit,saveplots=2,savedir);
screeplotfunc(datob=cctfit$datob,saveplots=1,savedir); screeplotfunc(datob=cctfit$datob,saveplots=2,savedir)
message("\n ...Export complete\n")
}
######################
#Function for the 'Fit Summary' Button
#- Gives a basic summary of the fit
######################
summaryfuncbutton <- function() {
guidat <- get("guidat", pkg_globals)
summary(get(guidat$varname, pkg_globals))
}
######################
#Function for the 'Print Fit' Button
#- Prints the cctfit object in the console
######################
printfitfuncbutton <- function() {
guidat <- get("guidat", pkg_globals)
print(get(guidat$varname, pkg_globals))
}
######################
#Function for the 'Send to Console' Button
#- Gives user the code to load the cctfit object into the console
######################
sendconsolefuncbutton <- function() {
guidat <- get("guidat", pkg_globals)
message("\n ...Paste the following to send your fit object into the R console
\n ",guidat$varname," <- CCTpack:::pkg_globals$",guidat$varname,
"\n",sep="")
}
######################
#Function for the 'Send to Console' Button
#- Gives user the code to load the cctfit object into the console
######################
sendconsolefuncbutton <- function() {
guidat <- get("guidat", pkg_globals)
message("",sep="")
message("\n ...Paste the following to send your fit object into the R console
\n ",guidat$varname," <- CCTpack:::pkg_globals$",guidat$varname,
"\n")
}
######################
#Function for the 'Resp Memberships' Button
#- Outputs the respondent memberships from the model fit
######################
membfuncbutton <- function(){
guidat <- get("guidat", pkg_globals)
message("\n ...Model-based cluster memberships for each respondent \n",sep="")
print(get(guidat$varname, pkg_globals)$respmem)
#print(get(guidat$varname, pkg_globals)$MVest)
}
######################
#Function for the 'NA Value Est' Button
#- Plot traceplots for all discrete parameters 3x3 plot size
######################
mvestfuncbutton <- function(){
guidat <- get("guidat", pkg_globals)
cctfit <- get(guidat$varname, pkg_globals)
if(cctfit$mval==FALSE){
message("\n ...Data has no missing values to estimate\n",sep="")
return()
}
if(cctfit$checksrun==FALSE){
message("\n ...Please use 'Run Checks' first",sep="")
}else{
message("\n ...Model estimates for missing values \n",sep="")
print(cctfit$MVest)
}
}
######################
#Function for the 'Traceplot Discrete' Button
#- Plot traceplots for all discrete parameters 3x3 plot size
######################
traceplotdiscretefuncbutton <- function() {
guidat <- get("guidat", pkg_globals)
dtraceplot(get(guidat$varname, pkg_globals))
}
######################
#Function for the 'Traceplot Continuous' Button
#- Plot traceplots for all parameters 3x3 plot size
######################
traceplotallfuncbutton <- function() {
guidat <- get("guidat", pkg_globals)
traceplot(get(guidat$varname, pkg_globals),mfrow=c(3,3), ask = FALSE)
}
dtraceplot <- function(cctfit,ask=FALSE){
if(cctfit$V==1 && cctfit$whmodel!="GCM"){
message("\n ...There are no discrete nodes in this inference")
return()
}
if(cctfit$whmodel!="GCM"){
inds <- c(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]])
}else{
if(cctfit$V==1){
inds <- c(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Z")]])
}
if(cctfit$V>1){
inds <- c(cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Z")]],cctfit$BUGSoutput$long.short[[which(cctfit$BUGSoutput$root.short=="Om")]])
}
}
if(cctfit$BUGSoutput$n.chains > 1){
cctfit$BUGSoutput$sims.array <- cctfit$BUGSoutput$sims.array[,,inds]
}else{
dmnames <- dimnames(cctfit$BUGSoutput$sims.array); dmnames[[3]] <- dmnames[[3]][inds]
tmp <- array(cctfit$BUGSoutput$sims.array[,,inds],c(cctfit$BUGSoutput$n.keep,length(inds),1))
cctfit$BUGSoutput$sims.array <- aperm(tmp,c(1,3,2)); dimnames(cctfit$BUGSoutput$sims.array) <- dmnames;
}
traceplot(cctfit,ask=ask,mfrow=c(3,3))
}
#######################
#End of Functions for the GUI and/or Standalone Functions
#######################
#######################
#Package Methods
#######################
setClass("cct",contains=c("rjags"));
setGeneric("plot"); setGeneric("screeplot"); setGeneric("summary")
setMethod("plot",signature=c(x="cct"),definition=function(x) cctresults(x))
setMethod("summary",signature=c(object="cct"),definition=function(object) cctsum(object))
setMethod("screeplot",signature=c(x="data.frame"),definition=function(x,polych=F) cctscree(x,polych=polych))
setMethod("screeplot",signature=c(x="matrix"),definition=function(x,polych=F) cctscree(x,polych=polych))
setMethod("screeplot",signature=c(x="cct"),definition=function(x,polych=F) cctscree(x$datob$dat,polych=polych))
######################
#END CCTpack
######################
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.