Bayes Analysis

Specify the DATA.

        y = CP_b[,"Pain.obsv"]
        x = CP_b[,"Pain.exp"]
        # Convert sName to consecutive integers:
        s = as.numeric(factor(CP_b[,"ID"]))
        # 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.") }
        #Ntotal = length(y)
        # Specify the data in a list, for later shipment to JAGS:
        dataList = list(
        x = x ,
        y = y ,
        s = s ,
        Nsubj = max(s)  # should equal length(unique(s))
        )

Specify the MODEL.

modelString = "
          # Standardize the data:
          data {
            Ntotal <- length(y)
            xm <- mean(x)
            ym <- mean(y)
            xsd <- sd(x)
            ysd <- sd(y)
            for ( i in 1:length(y) ) {
              zx[i] <- ( x[i] - xm ) / xsd
              zy[i] <- ( y[i] - ym ) / ysd
            }
          }
          # Specify the model for standardized data:
          model {
            for ( i in 1:Ntotal ) {
              zy[i] ~ dt( zbeta0[s[i]] + zbeta1[s[i]] * zx[i] , 1/zsigma^2 , nu )
            }
            for ( j in 1:Nsubj ) {
              zbeta0[j] ~ dnorm( zbeta0mu , 1/(zbeta0sigma)^2 )  
              zbeta1[j] ~ dnorm( zbeta1mu , 1/(zbeta1sigma)^2 )
            }
            # Priors vague on standardized scale:
            zbeta0mu ~ dnorm( 0 , 1/(10)^2 )
            zbeta1mu ~ dnorm( 0 , 1/(10)^2 )
            zsigma ~ dunif( 1.0E-3 , 1.0E+3 )
            zbeta0sigma ~ dunif( 1.0E-3 , 1.0E+3 )
            zbeta1sigma ~ dunif( 1.0E-3 , 1.0E+3 )
            nu ~ dexp(1/30.0)
            # Transform to original scale:
            for ( j in 1:Nsubj ) {
              beta1[j] <- zbeta1[j] * ysd / xsd  
              beta0[j] <- zbeta0[j] * ysd  + ym - zbeta1[j] * xm * ysd / xsd 
            }
            beta1mu <- zbeta1mu * ysd / xsd  
            beta0mu <- zbeta0mu * ysd  + ym - zbeta1mu * xm * ysd / xsd 
            sigma <- zsigma * ysd
          }
          " # close quote for modelString
          # Write out modelString to a text file
          writeLines( modelString , con="TEMPmodel.txt" ) #Write the model to a text file

RUN THE CHAINS

          numSavedSteps=10000 
          thinSteps = 1 
          saveName="JAGS_OUT"
          runjagsMethod=runjagsMethodDefault 
          nChains=3

          parameters = c( "beta0" ,  "beta1" , "beta0mu" , "beta1mu" ,
                          "zbeta0" , "zbeta1" ,  "zbeta0mu" , "zbeta1mu" ,
                          "zsigma", "sigma", "nu" , "zbeta0sigma" , "zbeta1sigma" )
          adaptSteps = 1000  # Number of steps to "tune" the samplers
          burnInSteps = 2000
          runJagsOut <- run.jags( method=runjagsMethod ,
                                  model="TEMPmodel.txt" , 
                                  monitor=parameters , 
                                  data=dataList ,  
                                  #inits=initsList , 
                                  n.chains=nChains ,
                                  adapt=adaptSteps ,
                                  burnin=burnInSteps , 
                                  sample=ceiling(numSavedSteps/nChains) ,
                                  thin=thinSteps ,
                                  summarise=FALSE ,
                                  plots=FALSE )
          codaSamples = as.mcmc.list( runJagsOut )
          # resulting codaSamples object has these indices: 
          #   codaSamples[[ chainIdx ]][ stepIdx , paramIdx ]
          if ( !is.null(saveName) ) {
            save( codaSamples , file=paste(saveName,"Mcmc.Rdata",sep="") )
          }

Load a Saved MCMC Chain

load("JAGS_OUTMcmc.Rdata")
mcmcCoda <- codaSamples

Diagnostics

# Diagonostics
windows.options(reset = TRUE)
parameterNames = varnames(mcmcCoda) # get all parameter names
for ( parName in c("beta0mu","beta1mu","nu","sigma","beta0[1]","beta1[1]") ) {
 diag.Plot <- diagMCMC( codaObject=mcmcCoda , parName=parName , 
            saveName=NULL , saveType=NULL )
 diag.Plot
}

Summary Info

summaryInfo = smryMCMC( mcmcCoda , saveName=NULL )
#show(summaryInfo)

windows.options(reset = TRUE)
plotMCMC( mcmcCoda , data=CP_b , xName="Pain.exp" , yName="Pain.obsv" , sName="ID" ,
          compValBeta1=0.0 , ropeBeta1=c(-0.5,0.5) ,
          pairsPlot=TRUE , showCurve=FALSE ,
          saveName=NULL , saveType=NULL )