library("rjags") source("coda2.R") #converts JAGS output to coda source("sortmix.R") #Has initial functions and #output sorting code source("mixfuns.R") #calculates WAIC and DIC ## Read in raw data source("simdata.R") ## Script file for running JAGS ## This selects the data file DataSet <- K4simdat DataSetName <- "K4 Simulation JAGS" ## This selects the essays whichLevel2 <- 1:10 ## Minimum number of entries before we will include it. This should ## be at least 1, and long enough to filter out "null" essays. minLen <- 25 ## Number of components to fit K <- 4 runName <- make.names(paste(DataSetName,K)) ## Control of plotting display plotParams <- TRUE printParams <- TRUE saveGraphs <- TRUE ## Parameters for MCMC n.chains <- 3 chain.frac <- c(1,.8,.8) # Chain starting points will # be calculated from a sample # of this fraction of each # data vector. chain.boot <- FALSE # Should this be a bootstrap # sample (drawn with # replacement)? n.retry <- 5 # Number of times to retry the # sampling before giving up. n.adapt <- 2000 n.burnin <- 1000 ## after adaptation. n.iter.level2 <- 5000 n.iter.level1 <- 2000 n.thin.level2 <- 1 n.thin.level1 <- 2 secondChance <- 2 #Length multiplier in case #first run does not converge mcmc.maxruns <- 2 #maximum number of MCMC runs to try ## Convergence criteria Rhat.max <- 1.2 neff.min <- 100 ## Parameters to monitor. Might need a subset of these. ## per component hyperparameters level2.params <- c("alpha0","mu0","tau0","beta0","gamma0") level2.scalarparams <- c("deviance","alphaN") # These parameters don't # depend on K, so don't need sorting. ## parameter which will be used to sort mixture components ## mu0 and alpha are sensible choices. May want to set decreasing to ## true with alpha. jags.model <- "hierModel1p.jags" level2.idparam <- "mu0" decreasing <- FALSE ## per essay/component parameteters level1.params <- c("pi","mu","tau") level1.idparam <- "mu" # This should correspond to # the level 2 parameter simplex.params <- c("pi","alpha0") #These params are defined on a #simplex, so requires some #care when using gelman.diag #function ## Level 2 hyperparameters mu0m <- 0 #location of mean log group location prior mu0s <- 1000 #scale of mean log group location prior lbeta0m <- 0 #location of sd log group location prior lbeta0s <- 1 #location of sd log group location prior tau0m <- 0 #location of mean log group precision prior tau0s <- 1 #scale of mean log group precision prior lgamma0m <- 0 #location of sd log group precision prior lgamma0s <- 1 #location of sd log group precision prior ## Set details of model in result object. This will be saved at the ## end of the day. result1 <- list(dataSetName=DataSetName,K=K,model=jags.model, n.iter.level1=n.iter.level1,n.adapt=n.adapt, n.iter.level2=n.iter.level2,n.burnin=n.burnin, n.thin.level2=n.thin.level2,n.thin.level1=n.thin.level1, constrained=FALSE, #Can't do this in JAGS idparam=level2.idparam,decreasing=decreasing) ## Wrap entire script in a giant tryCatch. The finally clause should ## cause a partial result1 to get written out for debugging in the ## case of a crash somewhere along the line tryCatch({ ### Start of actual script cat("\n\n ****************\n") cat("Cleaning data for ",DataSetName,"@",K,"\n") ## Extract desired data set rawdata<-list() rawdata$Y <- DataSet$Y[whichLevel2,] rawdata$J <- DataSet$J[whichLevel2] ## clean excessively short logs short <- rawdata$J < minLen cat("Removing ", sum(short), "of", length(short),"essays for length.\n") result1$data <- list(Y=rawdata$Y[!short,], J=rawdata$J[!short],I=sum(!short), alphaNu=sum(!short),K=K, alpha0mu=rep(0.5,K), mu0m=rep(mu0m,K), mu0s=rep(mu0s,K), lbeta0m=rep(lbeta0m,K), lbeta0s=rep(lbeta0s,K), tau0m=rep(tau0m,K), tau0s=rep(tau0s,K), lgamma0m=rep(lgamma0m,K), lgamma0s=rep(lgamma0s,K) ) ##Jmax is never actually used in the JAGS code ##result1$data$Jmax <- max(result1$data$J) #Precalc for array allocation ### Calculate Initial Values result1$inits <- vector("list",n.chains) for (d in 1:n.chains) { cat("Calculating initial values for chain ",d,";", DataSetName,"@",K,"\n") result1$inits[[d]] <- hmmInits(d,result1$data$Y,result1$data$J, n.retry=n.retry, decreasing=decreasing) } ## We could pass the hmmInits funciton instead of result1$inits, but this ## method allows us to save the inital values we used. ###Model Setup, Adapatation and Burnin cat("Loading Model for ",DataSetName,"@",K,"\n") load.module("mix") load.module("dic") heirModel <- jags.model(jags.model,data=result1$data, inits=result1$inits, n.chains=n.chains,n.adapt=n.adapt) ## Burn in cat("Burn in iterations for ",DataSetName,"@",K,"\n") update(heirModel,n.burnin) ### Hyperparameter learning loop. ## Using JAGS we will do a relatively long run monitoring the ## hyperparameters, and then a shorter run monitoring the parmaeters ## This code is probably redundant, but we will use this to test ## later for the first successful run result1$hparams <- NULL result1$scalars <- NULL result1$params <- NULL result1$failedruns <- numeric() result1$converged <- logical() ## We will put this into a while loop which will allow us to ## continue to the second chance if the first one fails to start for ## some reason converged <- FALSE runno <- 0 while (!converged && runno < mcmc.maxruns) { runno <- runno+1 if (runno > 1) { ## Double run length second time through n.iter.level2 <- n.iter.level2*secondChance result1$n.iter.level2 <-c(result1$n.iter.level2,n.iter.level2) } ## Monitor hyperparameters cat("\n\n\n ****************\n") cat("Learning hyperparameters for ",DataSetName,"@",K, "Attempt", runno,"\n") run1.samples <- jags.samples(heirModel,c(level2.params,level2.scalarparams), n.iter.level2,thin=n.thin.level2) ## Extract Level 2 and deviance parameters. run1.level2 <- jags2coda2(run1.samples,level2.params) run1.scalars <- jags2coda2(run1.samples,level2.scalarparams) ## Stan runs occasionally crash in the middle, returning ## essentially a null fit object. In this case, the extracted ## values should all be null. I haven't seen JAGS crash like ## this, but it might. if (is.null(run1.level2[[1]])) { result1$failedruns <- c(result1$failedruns,runno) cat("Run ",runno," crashed; skipping convergence checks.\n") next } ### Identify components ## The dmixnorm() function in JAGS does efficient sampling in the ## mixture distribution, which is multimodal, but does not identify ## the compoennts. The easiest way to identify the components is to ## pick a variable, e.g., mu0 or alpha, and identify component 1 as ## the mode with the smallest value of the variable, component 2 as ## the second smallest and so forth. cat("Labeling components for level 2 model",DataSetName,"@",K,"\n") ## Pull into MCMC list object, as MCarray object is broken. At the ## same time reorder the components into the order specified above. run1.level2 <- level2sort(run1.level2,level2.params, sortindex(run1.level2[[level2.idparam]], decreasing)) ## bind the sorted components onto the results if (is.null(result1$hparams)) { result1$hparams <- run1.level2 } else { for (param in level2.params) { result1$hparams[[param]] <- mcmcjoin(result1$hparams[[param]], run1.level2[[param]]) } } if (is.null(result1$scalars)) { result1$scalars <- run1.scalars } else { for (param in level2.scalarparams) { result1$scalars[[param]] <- mcmcjoin(result1$scalars[[param]], run1.scalars[[param]]) } } ### Convergence Diagnostics and printout cat("\n\n ****************\n") cat("Convergence diagnostics for ",DataSetName,"@",K, "Run Number", runno,"\n") ## Accumulate these as we go along, to test for convergence Rhatmax <- 0 neffMin <- Inf ## If we have null parameters, that means something went wrong with ## MCMC sampler, we will want to try again anyway. tryCatch({ if (saveGraphs) pdf(paste(runName,runno,"pdf",sep=".")) for (param in level2.scalarparams) { if (printParams) { print(try(summary(result1$scalars[[param]]))) Rhat <-try(gelman.diag(result1$scalars[[param]],autoburnin=FALSE)) print(Rhat) if (!inherits(Rhat, "try-error")) Rhatmax <- max(Rhatmax,Rhat$psrf[,1]) neff <- try(effectiveSize(result1$scalars[[param]])) print(neff) if (!inherits(Rhat, "try-error")) neffMin<- min(neffMin,neff) } } if (plotParams) { ## set up looping so gelman plots and trace plots appear on ## the same page. par(mfrow=c(length(level2.scalarparams),1)) for (param in level2.scalarparams) { try(gelman.plot(result1$scalars[[param]],auto.layout=FALSE)) } par(mfrow=c(length(level2.scalarparams),2)) for (param in level2.scalarparams) { try(plot(result1$scalars[[param]],auto.layout=FALSE)) } } for (param in level2.params) { if (printParams) { print(try(summary(result1$hparams[[param]]))) if (param %in% simplex.params) { ## Drop first column to avoid divide by zero error Rhat <-try(gelman.diag(result1$hparams[[param]][,-1], autoburnin=FALSE)) } else { Rhat <-try(gelman.diag(result1$hparams[[param]],autoburnin=FALSE)) } print(Rhat) if (!inherits(Rhat, "try-error")) Rhatmax <- max(Rhatmax,Rhat$mpsrf) neff <- try(effectiveSize(result1$hparams[[param]])) print(neff) if (!inherits(Rhat, "try-error")) neffMin<- min(neffMin,neff) } if (plotParams) { par(mfrow=c(K,1)) if (param %in% simplex.params) { ## Drop first column to avoid divide by zero error try(gelman.plot(result1$hparams[[param]][,-1],auto.layout=FALSE)) } else { try(gelman.plot(result1$hparams[[param]],auto.layout=FALSE)) } par(mfrow=c(K,1)) try(plot(result1$hparams[[param]])) } } }, finally=if (saveGraphs) dev.off()) ### Automatic convergence check converged <- TRUE if (Rhatmax > Rhat.max) { converged <- FALSE cat("Chains of length ",n.iter.level2,"for ", DataSetName,"@",K,"did not converge in run ",runno,".\n") cat("Maximum Rhat value = ",Rhatmax,".\n") if (printParams) { ## Print some hopefully diagnostic information about what the ## chains look like for (param in level2.scalarparams) { for (chain in 1:length(result1$scalars[[param]])) { cat(param,"[[",chain,"]]\n") print(try(summary(result1$scalars[[param]][[chain]])$stat)) } } for (param in level2.params) { for (chain in 1:length(result1$hparams[[param]])) { cat(param,"[[",chain,"]]\n") print(try(summary(result1$hparams[[param]][[chain]])$stat)) } } } } if (converged && neffMin < neff.min) { cat("Chain length of ",n.iter.level2,"for ", DataSetName,"@",K,"were insufficient in run",runno,".\n") cat("Minimum effective MC sample size = ",neffMin,".\n") converged <- FALSE } result1$converged <- c(result1$converged,converged) cat("\n\n\n") } ### End of Run loop. if (!converged) { cat("MCMC run did not converge, proceeding anyway.\n") } ### Additional 3rd run for parameters and DIC ## Monitor parameters cat("Learning parameters for ",DataSetName,"@",K,"\n") run3.samples <- jags.samples(heirModel,c(level2.idparam,level1.params), n.iter.level1,thin=n.thin.level1) ## Extract parameters run3.idparam <- jags2coda(run3.samples,level2.idparam) if (is.null(run3.idparam)) { stop("Level 1 run failed for some reason.") } run3.level1 <- jags2coda2(run3.samples,level1.params,K) ### Identify components on second level parameters ## Monitor component assignments Need to do component assignments ## post hoc. cat("Labeling components for ",DataSetName,"@",K,"\n") run3.index <- sortindex(run3.idparam,decreasing) result1$params <- level1sort(run3.level1,level1.params,run3.index,K) cat("Calculating model fit indexes for ",DataSetName,"@",K,"\n") sample.pi <- Karray(combineChains(result1$params$pi),K) sample.mu <- Karray(combineChains(result1$params$mu),K) sample.sigma <- 1/sqrt(Karray(combineChains(result1$params$tau),K)) result1$WAIC <- mchmixWAIC(result1$data$Y,result1$data$J, sample.pi,sample.mu,sample.sigma) print(result1$WAIC) result1$DIC <- mchmixDIC(result1$data$Y,result1$data$J, sample.pi,sample.mu,sample.sigma) print(result1$DIC) cat("Analaysis complete for ",DataSetName,"@",K,"\n") }, error=function(e) print(e), finally= { ### Write out results ## We do this in the finally clause because it that way the ## partial object will get written even if there is an ## error above. assign(runName,result1) outfile <- gzfile(paste(runName,"R","gz",sep="."),open="wt") dump(runName,file=outfile) close(outfile) } )