inst/examples/Stan-Ydich-XnomSsubj-MbernBetaOmegaKappa-Example.R

In variani/dbda: R code from Doing Bayesian Data Analysis (DBDA)

# Example for Stan-Ydich-XnomSsubj-MbernBetaOmegaKappa.R 
#------------------------------------------------------------------------------- 
# Optional generic preliminaries:
#graphics.off() # This closes all of R's graphics windows.
#rm(list=ls())  # Careful! This clears all of R's memory!
#------------------------------------------------------------------------------- 
# Load The data 
#myData = read.csv("TherapeuticTouchData.csv")
myData = read.csv(system.file("data", "TherapeuticTouchData.csv", package = "dbda"))
# N.B.: The functions below expect 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.
#------------------------------------------------------------------------------- 
# Load the relevant model into R's working memory:
#source("Stan-Ydich-XnomSsubj-MbernBetaOmegaKappa.R")
#------------------------------------------------------------------------------- 
# Optional: Specify filename root and graphical format for saving output.
# Otherwise specify as NULL or leave saveName and saveType arguments 
# out of function calls.
#fileNameRoot = "Stan-Ydich-XnomSsubj-MbernBetaOmegaKappa-" 
fileNameRoot = NULL
#graphFileType = "eps" 
#------------------------------------------------------------------------------- 


StanYdichXnomSsubjMbernBetaOmegaKappa <- function()
{
  out <- list()
  
  oldClass(out) <- "StanYdichXnomSsubjMbernBetaOmegaKappa"
  return(out)
}

mod <- StanYdichXnomSsubjMbernBetaOmegaKappa()

# Generate the MCMC chain:
startTime = proc.time()
out = genMCMC(mod, data=myData , sName="s" , yName="y" ,  
                    numSavedSteps=12000 , saveName=fileNameRoot , thinSteps=10 )

mcmcCoda <- out$codaSamples
stanFit <- out$stanFit           

stopTime = proc.time()
duration = stopTime - startTime
show(duration)

#------------------------------------------------------------------------------- 
# Display diagnostics of chain, for specified parameters:
parameterNames = varnames(mcmcCoda) # get all parameter names for reference
for ( parName in c("omega","kappa","theta[1]","theta[2]") ) { 
  diagMCMC( codaObject=mcmcCoda , parName=parName , 
                saveName=fileNameRoot , saveType=graphFileType )
}
#------------------------------------------------------------------------------- 
# Get summary statistics of chain:
summaryInfo = smryMCMC(mod, mcmcCoda , compVal=0.5 , 
                        diffIdVec=c(1,14,28), compValDiff=0.0, 
                        saveName=fileNameRoot )
# Display posterior information:
plotMCMC(mod, mcmcCoda , data=myData , compVal=0.5 , #rope=c(0.45,0.55) ,
          diffIdVec=c(1,14,28), compValDiff=0.0, #ropeDiff = c(-0.05,0.05) ,
          saveName=fileNameRoot , saveType=graphFileType )
#------------------------------------------------------------------------------- 
# Using Stan display functions instead of DBDA2E functions:
# Load the stanFit object that was saved by genMCMC:
#load(paste0(fileNameRoot,"StanFit.Rdata"))
# Display information:
show(stanFit)
#openGraph()
traceplot(stanFit,pars= c("omega","kappa","theta[1]","theta[2]"))
#openGraph()
plot(stanFit)