scripts update

2e/DBDA2E-utilities.R

@@ -1,38 +1,35 @@
 # Utility programs for use with the book,
-#   Kruschke, J. K. (2014). Doing Bayesian Data Analysis: 
-#   A Tutorial with R, JAGS, and Stan. 2nd Edition. Academic Press.
+#   Kruschke, J. K. (2015). Doing Bayesian Data Analysis, Second Edition: 
+#   A Tutorial with R, JAGS, and Stan. Academic Press / Elsevier.
 # This file contains several functions that are called by other programs
 # or can be called directly by the user. To load all the functions into
 # R's working memory, at R's command line type:
 # source("DBDA2E-utilities.R")
 
 #------------------------------------------------------------------------------
 
-bookInfo = "Kruschke, J. K. (2014). Doing Bayesian Data Analysis:\nA Tutorial with R, JAGS, and Stan. 2nd Edition. Academic Press."
-bannerBreak = "\n*******************************************************************\n"
+bookInfo = "Kruschke, J. K. (2015). Doing Bayesian Data Analysis, Second Edition:\nA Tutorial with R, JAGS, and Stan. Academic Press / Elsevier."
+bannerBreak = "\n*********************************************************************\n"
 cat(paste0(bannerBreak,bookInfo,bannerBreak,"\n"))
 
 #------------------------------------------------------------------------------
 # Check that required packages are installed:
+want = c("parallel","rjags","runjags")
+have = want %in% rownames(installed.packages())
+if ( any(!have) ) { install.packages( want[!have] ) }
 
-pkgTest <- function( pkgName ) {
-  if ( !require( pkgName , character.only=TRUE ) ) {
-    install.packages( pkgName , dep=TRUE )
-    if ( !require( pkgName , character.only=TRUE ) ) { 
-      stop(paste0("Package ", pkgName ," not found")) 
-    }
-  }
-}
-
-requiredPackages = c("parallel","rjags","runjags") #  rstan is not on CRAN!
-for ( pkgName in requiredPackages ) { pkgTest( pkgName ) }
+# Load parallel package for detectCores function:
+library(parallel)
+# Set some options in runjags:
+library(runjags)
+runjags.options( inits.warning=FALSE , rng.warning=FALSE )
 
 #------------------------------------------------------------------------------
 
 # Make some data files for examples...
 createDataFiles=FALSE
 if ( createDataFiles ) {
-
+  
   source("HtWtDataGenerator.R")
   N=300
   m = HtWtDataGenerator( N , rndsd=47405 )
@@ -78,13 +75,24 @@ if ( createDataFiles ) {
 
 #------------------------------------------------------------------------------
 # Functions for opening and saving graphics that operate the same for 
-# Windows and Macintosh and Linux operating systems. At least, that's the hope.
+# Windows and Macintosh and Linux operating systems. At least, that's the hope!
 
 openGraph = function( width=7 , height=7 , mag=1.0 , ... ) {
   if ( .Platform$OS.type != "windows" ) { # Mac OS, Linux
-    X11( width=width*mag , height=height*mag , type="cairo" , ... )
+    tryInfo = try( X11( width=width*mag , height=height*mag , type="cairo" , 
+                        ... ) )
+    if ( class(tryInfo)=="try-error" ) {
+      lineInput = readline("WARNING: Previous graphics windows will be closed because of too many open windows.\nTO CONTINUE, PRESS <ENTER> IN R CONSOLE.\n")
+      graphics.off() 
+      X11( width=width*mag , height=height*mag , type="cairo" , ... )
+    }
   } else { # Windows OS
-    windows( width=width*mag , height=height*mag , ... )
+    tryInfo = try( windows( width=width*mag , height=height*mag , ... ) )
+    if ( class(tryInfo)=="try-error" ) {
+      lineInput = readline("WARNING: Previous graphics windows will be closed because of too many open windows.\nTO CONTINUE, PRESS <ENTER> IN R CONSOLE.\n")
+      graphics.off() 
+      windows( width=width*mag , height=height*mag , ... )    
+    }
   }
 }
 
@@ -233,7 +241,7 @@ DbdaDensPlot = function( codaObject , parName=varnames(codaObject)[1] , plColors
 }
 
 diagMCMC = function( codaObject , parName=varnames(codaObject)[1] ,
-                         saveName=NULL , saveType="jpg" ) {
+                     saveName=NULL , saveType="jpg" ) {
   DBDAplColors = c("skyblue","black","royalblue","steelblue")
   openGraph(height=5,width=7)
   par( mar=0.5+c(3,4,1,0) , oma=0.1+c(0,0,2,0) , mgp=c(2.25,0.7,0) , 
@@ -246,12 +254,17 @@ diagMCMC = function( codaObject , parName=varnames(codaObject)[1] ,
   tryVal = try(
     coda::gelman.plot( codaObject[,c(parName)] , main="" , auto.layout=FALSE , 
                        col=DBDAplColors )
-  )
+  )  
   # if it runs, gelman.plot returns a list with finite shrink values:
-  if ( class(tryVal) != "list"  | !is.finite(tryVal$shrink[1]) ) { 
+  if ( class(tryVal)=="try-error" ) {
     plot.new() 
-    print(paste0("Warning: gelman.plot fails for ",parName))
-  } 
+    print(paste0("Warning: coda::gelman.plot fails for ",parName))
+  } else { 
+    if ( class(tryVal)=="list" & !is.finite(tryVal$shrink[1]) ) {
+      plot.new() 
+      print(paste0("Warning: coda::gelman.plot fails for ",parName))
+    }
+  }
   DbdaAcfPlot(codaObject,parName,plColors=DBDAplColors)
   DbdaDensPlot(codaObject,parName,plColors=DBDAplColors)
   mtext( text=parName , outer=TRUE , adj=c(0.5,0.5) , cex=2.0 )
@@ -261,7 +274,7 @@ diagMCMC = function( codaObject , parName=varnames(codaObject)[1] ,
 }
 
 diagStanFit = function( stanFit , parName ,
-                     saveName=NULL , saveType="jpg" ) {
+                        saveName=NULL , saveType="jpg" ) {
   codaFit = mcmc.list( lapply( 1:ncol(stanFit) , 
                                function(x) { mcmc(as.array(stanFit)[,x,]) } ) )
   DBDAplColors = c("skyblue","black","royalblue","steelblue")
@@ -274,9 +287,19 @@ diagStanFit = function( stanFit , parName ,
   # gelman.plot are from CODA package:
   require(coda)
   tryVal = try(
-    gelman.plot(codaFit[,c(parName)],main="",auto.layout=FALSE,col=DBDAplColors)
+    coda::gelman.plot( codaObject[,c(parName)] , main="" , auto.layout=FALSE , 
+                       col=DBDAplColors )
   )
-  if ( class(tryVal)!="list" ) { plot.new() } # gelman.plot returns list
+  # if it runs, gelman.plot returns a list with finite shrink values:
+  if ( class(tryVal)=="try-error" ) {
+    plot.new() 
+    print(paste0("Warning: coda::gelman.plot fails for ",parName))
+  } else { 
+    if ( class(tryVal)=="list" & !is.finite(tryVal$shrink[1]) ) {
+      plot.new() 
+      print(paste0("Warning: coda::gelman.plot fails for ",parName))
+    }
+  }
   DbdaAcfPlot(codaFit,parName,plColors=DBDAplColors)
   DbdaDensPlot(codaFit,parName,plColors=DBDAplColors)
   mtext( text=parName , outer=TRUE , adj=c(0.5,0.5) , cex=2.0 )
@@ -427,16 +450,16 @@ plotPost = function( paramSampleVec , cenTend=c("mode","median","mean")[1] ,
            lty="dotted" , lwd=2 , col=cvCol )
     text( compVal , cvHt ,
           bquote( .(round(100*pLtCompVal,1)) * "% < " *
-                  .(signif(compVal,3)) * " < " * 
-                  .(round(100*pGtCompVal,1)) * "%" ) ,
+                   .(signif(compVal,3)) * " < " * 
+                   .(round(100*pGtCompVal,1)) * "%" ) ,
           adj=c(pLtCompVal,0) , cex=0.8*cex , col=cvCol )
     postSummary[,"compVal"] = compVal
     postSummary[,"pGtCompVal"] = pGtCompVal
   }
   # Display the ROPE.
   if ( !is.null( ROPE ) ) {
     pInROPE = ( sum( paramSampleVec > ROPE[1] & paramSampleVec < ROPE[2] )
-                 / length( paramSampleVec ) )
+                / length( paramSampleVec ) )
     pGtROPE = ( sum( paramSampleVec >= ROPE[2] ) / length( paramSampleVec ) )
     pLtROPE = ( sum( paramSampleVec <= ROPE[1] ) / length( paramSampleVec ) )
     lines( c(ROPE[1],ROPE[1]) , c(0.96*ROPEtextHt,0) , lty="dotted" , lwd=2 ,
@@ -445,8 +468,8 @@ plotPost = function( paramSampleVec , cenTend=c("mode","median","mean")[1] ,
            col=ropeCol)
     text( mean(ROPE) , ROPEtextHt ,
           bquote( .(round(100*pLtROPE,1)) * "% < " * .(ROPE[1]) * " < " * 
-                  .(round(100*pInROPE,1)) * "% < " * .(ROPE[2]) * " < " * 
-                  .(round(100*pGtROPE,1)) * "%" ) ,
+                   .(round(100*pInROPE,1)) * "% < " * .(ROPE[2]) * " < " * 
+                   .(round(100*pGtROPE,1)) * "%" ) ,
           adj=c(pLtROPE+.5*pInROPE,0) , cex=1 , col=ropeCol )
 
     postSummary[,"ROPElow"]=ROPE[1] 

2e/Jags-Ymet-XmetSsubj-MrobustHier.R

@@ -15,7 +15,8 @@ genMCMC = function( data , xName="x" , yName="y" , sName="s" ,
   # THE DATA.
   y = data[,yName]
   x = data[,xName]
-  s = as.numeric(data[,sName])
+  # Convert sName to consecutive integers:
+  s = as.numeric(factor(data[,sName]))
   # Do some checking that data make sense:
   if ( any( !is.finite(y) ) ) { stop("All y values must be finite.") }
   if ( any( !is.finite(x) ) ) { stop("All x values must be finite.") }
@@ -242,9 +243,9 @@ plotMCMC = function( codaSamples , data , xName="x" , yName="y" , sName="s" ,
       }
       points( x[thisSrows] , y[thisSrows] , pch=19 )
     }
-  }
-  if ( !is.null(saveName) ) {
-    saveGraph( file=paste(saveName,"PostPredSubj",sep=""), type=saveType)
+    if ( !is.null(saveName) ) {
+      saveGraph( file=paste0(saveName,"PostPredSubj",plotIdx), type=saveType)
+    }
   }
   #-----------------------------------------------------------------------------
   # Data with superimposed regression lines and noise distributions:

2e/Jags-Ymet-XmetSsubj-MrobustHierQuadWt.R

@@ -103,7 +103,7 @@ genMCMC = function( data , xName="x" , yName="y" , sName="s" , wName=NULL ,
   nuInit = 10 # arbitrary
   initsList = list(
     zsigma=sigmaInit  ,
-    nu=nuInit ,
+    nuMinusOne=nuInit ,
     zbeta0mu=b0init ,
     zbeta1mu=b1init ,
     zbeta2mu=b2init ,

2e/Jags-Ymet-Xnom2fac-MnormalHom.R

@@ -343,7 +343,6 @@ plotMCMC = function( codaSamples ,
     if ( !is.null(saveName) ) {
       saveGraph( file=paste0(saveName,"PostPred-",x2levels[x2idx]), type=saveType)
     }
-    if ( length(x2levels) > 10 ) { dev.off() }
   } # end for x2idx
 
   # plot contrasts:

2e/Jags-Ymet-Xnom2fac-MrobustHet.R

@@ -359,7 +359,6 @@ plotMCMC = function( codaSamples ,
     if ( !is.null(saveName) ) {
       saveGraph( file=paste0(saveName,"PostPred-",x2levels[x2idx]), type=saveType)
     }
-    dev.off()
   } # end for x2idx
 
   # Display the normality and top-level variance parameters:

2e/Jags-Ymet-XnomSplitPlot-MnormalHom-Example.R

@@ -34,7 +34,7 @@ source("Jags-Ymet-XnomSplitPlot-MnormalHom.R")
 # Generate the MCMC chain:
 mcmcCoda = genMCMC( datFrm=myDataFrame , 
                     yName=yName , xBetweenName=xBetweenName , 
-                    xWithinName=xWithinName , xSubjName=xSubjectName ,
+                    xWithinName=xWithinName , xSubjectName=xSubjectName ,
                     numSavedSteps=12000 , thinSteps=20 , saveName=fileNameRoot )
 #------------------------------------------------------------------------------- 
 # Display diagnostics of chain, for specified parameters:

2e/Jags-Ymet-XnomSplitPlot-MnormalHom.R

@@ -3,7 +3,7 @@ source("DBDA2E-utilities.R")
 
 #===============================================================================
 
-genMCMC = function( datFrm , yName , xBetweenName , xWithinName , xSubjName ,
+genMCMC = function( datFrm , yName , xBetweenName , xWithinName , xSubjectName ,
                     numSavedSteps=50000 , thinSteps=1 , saveName=NULL ) { 
   #------------------------------------------------------------------------------
   # THE DATA.

2e/Jags-YmetBinned-Xnom1grp-MnormalInterval.R

@@ -1,7 +1,7 @@
 # Jags-YmetBinned-Xnom1grp-MnormalInterval.R 
 # Accompanies the book:
-#   Kruschke, J. K. (2014). Doing Bayesian Data Analysis: 
-#   A Tutorial with R and JAGS, 2nd Edition. Academic Press / Elsevier.
+#   Kruschke, J. K. (2014). Doing Bayesian Data Analysis, Second Edition: 
+#   A Tutorial with R, JAGS, and Stan. Academic Press / Elsevier.
 
 source("DBDA2E-utilities.R")
 
@@ -15,6 +15,7 @@ genMCMC = function( data ,  yName , yBinName , threshVecName ,
   y = data[,yName]
   ybin = data[,yBinName]
   threshMat = data[,threshVecName]
+  if ( is.null(dim(threshMat)) ) { threshMat = cbind(threshMat) }
   Ntotal = length(y)
   # Specify the data in a list, for later shipment to JAGS:
   dataList = list(
@@ -161,6 +162,7 @@ plotMCMC = function( codaSamples , data , yName , yBinName , threshVecName ,
   # Compute limits for plots of data with posterior pred. distributions
   y = data[,yName]
   threshMat = data[,threshVecName]
+  if ( is.null(dim(threshMat)) ) { threshMat = cbind(threshMat) }
   xLim = c( min(y,na.rm=TRUE)-0.1*(max(y,na.rm=TRUE)-min(y,na.rm=TRUE)) , 
             max(y,na.rm=TRUE)+0.1*(max(y,na.rm=TRUE)-min(y,na.rm=TRUE)) )
   xBreaks = seq( xLim[1] , xLim[2] , 

2e/Stan-Ydich-XnomSsubj-MbernBetaOmegaKappa-Example.R

@@ -20,8 +20,12 @@ fileNameRoot = "Stan-Ydich-XnomSsubj-MbernBetaOmegaKappa-"
 graphFileType = "eps" 
 #------------------------------------------------------------------------------- 
 # Generate the MCMC chain:
-mcmcCoda = genMCMC( data=myData , numSavedSteps=12000 , saveName=fileNameRoot ,
-                    thinSteps=10 )
+startTime = proc.time()
+mcmcCoda = genMCMC( data=myData , sName="s" , yName="y" ,  
+                    numSavedSteps=12000 , saveName=fileNameRoot , thinSteps=10 )
+stopTime = proc.time()
+duration = stopTime - startTime
+show(duration)
 #------------------------------------------------------------------------------- 
 # Display diagnostics of chain, for specified parameters:
 parameterNames = varnames(mcmcCoda) # get all parameter names for reference

2e/Stan-Ydich-XnomSsubj-MbernBetaOmegaKappa.R

@@ -5,16 +5,17 @@
 source("DBDA2E-utilities.R")
 #===============================================================================
 
-genMCMC = function( data , numSavedSteps=50000 , saveName=NULL , thinSteps=1 ) { 
+genMCMC = function( data , sName="s" , yName="y" ,  
+                    numSavedSteps=50000 , saveName=NULL , thinSteps=1 ) { 
   require(rjags)
   require(rstan)
   #-----------------------------------------------------------------------------
   # THE DATA.
   # N.B.: This function expects the data to be a data frame, 
   # with one component named y being a vector of integer 0,1 values,
   # and one component named s being a factor of subject identifiers.
-  y = as.numeric(data$y)
-  s = as.numeric(data$s) # converts character to consecutive integer levels
+  y = as.numeric(data[,yName])
+  s = as.numeric(data[,sName]) # ensures consecutive integer levels
   # Do some checking that data make sense:
   if ( any( y!=0 & y!=1 ) ) { stop("All y values must be 0 or 1.") }
   Ntotal = length(y)