Jags-Yord-Xnom1grp-Mnormal.R (11947B)
1 # Note: this code has been pulled over to antivax-attitudes.Rmd and subsequently 2 # modified. I'm keeping the file around for archival purposes. 3 4 # Jags-Yord-Xnom1grp-Mnormal.R 5 # Accompanies the book: 6 # Kruschke, J. K. (2015). Doing Bayesian Data Analysis, Second Edition: 7 # A Tutorial with R, JAGS, and Stan. Academic Press / Elsevier. 8 9 source("DBDA2E-utilities.R") 10 11 #=============================================================================== 12 13 genMCMC = function( datFrm, yName , x1Name, x2Name, x3Name, 14 numSavedSteps=50000 , thinSteps = 1, 15 saveName=NULL , 16 runjagsMethod=runjagsMethodDefault , 17 nChains=nChainsDefault ) { 18 #----------------------------------------------------------------------------- 19 # THE DATA. 20 x1 = as.numeric(as.factor(datFrm[[x1Name]])) 21 Nx1Lvl = max(x1) 22 23 x2 = as.numeric(as.factor(datFrm[[x2Name]])) 24 Nx2Lvl = max(x2) 25 26 x3 = as.numeric(as.factor(datFrm[[x3Name]])) 27 Nx3Lvl = max(x3) 28 29 y = as.numeric(datFrm[[yName]]) 30 # Do some checking that data make sense: 31 if ( any( y!=round(y) ) ) { stop("All y values must be integers (whole numbers).") } 32 if ( any( y < 1 ) ) { stop("All y values must be 1 or larger.") } 33 # COMPRESS OUT ANY EMPTY VALUES OF Y: 34 yOrig=y 35 y=as.numeric(factor(y,levels=names(table(y)))) 36 if ( any(y != yOrig) ) { 37 warning("*** WARNING: Y RE-CODED TO REMOVE EMPTY LEVELS ***") 38 } 39 Ntotal = length(y) 40 # For hyper-prior on deflections: 41 agammaShRa = unlist(gammaShRaFromModeSD(mode=sd(y)/2 , sd=2*sd(y))) 42 43 # Threshold 1 and nYlevels-1 are fixed; other thresholds are estimated. 44 # This allows all parameters to be interpretable on the response scale. 45 nYlevels = max(y) 46 thresh = matrix(data = NA, nrow = Nx1Lvl, ncol = nYlevels-1) 47 thresh[, 1] = 1 + 0.5 48 thresh[, nYlevels-1] = nYlevels-1 + 0.5 49 # Specify the data in a list, for later shipment to JAGS: 50 dataList = list( 51 x1 = x1, 52 Nx1Lvl = Nx1Lvl, 53 x2 = x2, 54 Nx2Lvl = Nx2Lvl, 55 x3 = x3, 56 Nx3Lvl = Nx3Lvl, 57 y = y, 58 NyLvl = nYlevels, 59 thresh = thresh, 60 Ntotal = Ntotal, 61 agammaShRa = agammaShRa 62 ) 63 #----------------------------------------------------------------------------- 64 # THE MODEL. 65 modelString = " 66 model { 67 for (i in 1:Ntotal) { 68 # Thresholded cummulative normal distribution 69 y[i] ~ dcat(pr[i,1:NyLvl]) 70 pr[i,1] <- pnorm(thresh[x1[i], 1], mu[i], 1/sigma[x1[i]]^2) 71 for (k in 2:(NyLvl-1)) { 72 pr[i,k] <- max(0, pnorm(thresh[x1[i], k] , mu[i] , 1/sigma[x1[i]]^2 ) - 73 pnorm(thresh[x1[i], k-1] , mu[i] , 1/sigma[x1[i]]^2 )) 74 } 75 pr[i,NyLvl] <- 1 - pnorm(thresh[x1[i], NyLvl-1] , mu[i] , 1/sigma[x1[i]]^2) 76 77 # mu ~ x1*x2*x3 78 mu[i] <- a0 + a1[x1[i]] + a2[x2[i]] + a3[x3[i]] + 79 a1a2[x1[i], x2[i]] + a1a3[x1[i], x3[i]] + a2a3[x2[i], x3[i]] + 80 a1a2a3[x1[i], x2[i], x3[i]] 81 } 82 83 a0 ~ dnorm((1+NyLvl)/2, 1/(NyLvl)^2) 84 85 for (j1 in 1:Nx1Lvl) { 86 # Constant sigma for beta1, we're treating all Qs as independent 87 a1[j1] ~ dnorm(0.0, 1/(NyLvl)^2) 88 89 # Sigma for normal CDF, unique for each x1. 90 sigma[j1] ~ dunif(NyLvl/1000, NyLvl*10) 91 92 # Threshold distributions. 1 and NyLvl-1 are fixed, not stochastic 93 for (k in 2:(NyLvl-2)) { 94 thresh[j1, k] ~ dnorm(k+0.5, 1/2^2) 95 } 96 } 97 98 # Constant sigma for beta2, the interventions are independent 99 for (j2 in 1:Nx2Lvl) { 100 a2[j2] ~ dnorm(0.0, 1/(NyLvl)^2) 101 } 102 103 # Constant sigma for beta3 (for now) 104 for (j3 in 1:Nx3Lvl) { 105 a3[j3] ~ dnorm(0.0, 1/(NyLvl)^2) 106 } 107 108 # Interaction terms also has homogenous variance (for now) 109 for (j1 in 1:Nx1Lvl) { 110 for (j2 in 1:Nx2Lvl) { 111 a1a2[j1, j2] ~ dnorm(0.0, 1/(NyLvl)^2) 112 } 113 } 114 for (j1 in 1:Nx1Lvl) { 115 for (j3 in 1:Nx3Lvl) { 116 a1a3[j1, j3] ~ dnorm(0.0, 1/(NyLvl)^2) 117 } 118 } 119 for (j2 in 1:Nx2Lvl) { 120 for (j3 in 1:Nx3Lvl) { 121 a2a3[j2, j3] ~ dnorm(0.0, 1/(NyLvl)^2) 122 } 123 } 124 for (j1 in 1:Nx1Lvl) { 125 for (j2 in 1:Nx2Lvl) { 126 for (j3 in 1:Nx3Lvl) { 127 a1a2a3[j1, j2, j3] ~ dnorm(0.0, 1/(NyLvl)^2) 128 } 129 } 130 } 131 132 # Compute cell means 133 for (j1 in 1:Nx1Lvl) { 134 for (j2 in 1:Nx2Lvl) { 135 for (j3 in 1:Nx3Lvl) { 136 m[j1, j2, j3] <- a0 + a1[j1] + a2[j2] + a3[j3] + 137 a1a2[j1, j2] + a1a3[j1, j3] + a2a3[j2, j3] + 138 a1a2a3[j1, j2, j3] 139 } 140 } 141 } 142 143 # Convert a0, a1[], a2[], &c. to sum-to-zero b0, b1[], b2[], &c. 144 b0 <- mean(m[1:Nx1Lvl, 1:Nx2Lvl, 1:Nx3Lvl]) 145 for (j1 in 1:Nx1Lvl) { 146 b1[j1] <- mean(m[j1, 1:Nx2Lvl, 1:Nx3Lvl]) - b0 147 } 148 for (j2 in 1:Nx2Lvl) { 149 b2[j2] <- mean(m[1:Nx1Lvl, j2, 1:Nx3Lvl]) - b0 150 } 151 for (j3 in 1:Nx3Lvl) { 152 b3[j3] <- mean(m[1:Nx1Lvl, 1:Nx2Lvl, j3]) - b0 153 } 154 for (j1 in 1:Nx1Lvl) { 155 for (j2 in 1:Nx2Lvl) { 156 b1b2[j1, j2] <- mean(m[j1, j2, 1:Nx3Lvl]) - (b0 + b1[j1] + b2[j2]) 157 } 158 } 159 for (j1 in 1:Nx1Lvl) { 160 for (j3 in 1:Nx3Lvl) { 161 b1b3[j1, j3] <- mean(m[j1, 1:Nx2Lvl, j3]) - (b0 + b1[j1] + b3[j3]) 162 } 163 } 164 for (j2 in 1:Nx2Lvl) { 165 for (j3 in 1:Nx3Lvl) { 166 b2b3[j2, j3] <- mean(m[1:Nx1Lvl, j2, j3]) - (b0 + b2[j2] + b3[j3]) 167 } 168 } 169 for (j1 in 1:Nx1Lvl) { 170 for (j2 in 1:Nx2Lvl) { 171 for (j3 in 1:Nx3Lvl) { 172 b1b2b3[j1, j2, j3] <- m[j1, j2, j3] - (b0 + b1[j1] + b2[j2] + b3[j3] + 173 b1b2[j1, j2] + b1b3[j1, j3] + b2b3[j2, j3]) 174 } 175 } 176 } 177 } 178 " # close quote for modelString 179 # Write out modelString to a text file 180 writeLines( modelString , con="TEMPmodel.txt" ) 181 #----------------------------------------------------------------------------- 182 # This is where the chains would be initialized, but we'll just let JAGS do it 183 initsList = NULL 184 #----------------------------------------------------------------------------- 185 # RUN THE CHAINS 186 parameters = c("b0", "b1", "b2", "b3", "b1b2", "b1b3", "b2b3", "b1b2b3", 187 "sigma", "thresh") 188 adaptSteps = 500 # Number of steps to "tune" the samplers 189 burnInSteps = 1000 190 runJagsOut <- run.jags( method=runjagsMethod , 191 model="TEMPmodel.txt" , 192 monitor=parameters , 193 data=dataList , 194 #inits=initsList , 195 n.chains=nChains , 196 adapt=adaptSteps , 197 burnin=burnInSteps , 198 sample=ceiling(numSavedSteps/nChains) , 199 thin=thinSteps , 200 summarise=FALSE , 201 plots=FALSE ) 202 codaSamples = as.mcmc.list( runJagsOut ) 203 # resulting codaSamples object has these indices: 204 # codaSamples[[ chainIdx ]][ stepIdx , paramIdx ] 205 if ( !is.null(saveName) ) { 206 save( codaSamples , file=paste(saveName,"Mcmc.Rdata",sep="") ) 207 } 208 return( codaSamples ) 209 } # end function 210 211 #=============================================================================== 212 213 plotMCMC = function( codaSamples , datFrm , yName , qName, compVal , #RopeEff=NULL , 214 showCurve=FALSE , 215 saveName=NULL , saveType="jpg" ) { 216 #----------------------------------------------------------------------------- 217 mcmcMat = as.matrix(codaSamples,chains=TRUE) 218 chainLength = NROW( mcmcMat ) 219 x1 = as.numeric(as.factor(datFrm[[qName]])) 220 Nx1Lvl = max(x1) 221 222 #----------------------------------------------------------------------------- 223 # Plots for each question 224 for (i in 1:Nx1Lvl) { 225 mu = mcmcMat[, "b0"] + mcmcMat[, paste0("b1[", i, "]")] 226 sigma = mcmcMat[, paste0("sigma[", i, "]")] 227 228 # Set up window and layout: 229 openGraph(width=6.0,height=5.0) 230 layout( matrix( c(1,2,3,4,5,6) , nrow=3, byrow=TRUE ) ) 231 par( mar=c(3.5,3.5,2.5,0.5) , mgp=c(2.25,0.7,0) ) 232 233 # Compute limits for plots of data with posterior pred. distributions 234 y = datFrm[[yName]][x1 == i] 235 xLim = c( min(y)-0.5 , max(y)+0.5 ) 236 xBreaks = seq( xLim[1] , xLim[2] , 1 ) 237 histInfo = hist(y,breaks=xBreaks,plot=FALSE) 238 yMax = 1.2 * max( c( histInfo$density ) ) 239 xVec = seq( xLim[1] , xLim[2] , length=501 ) 240 241 #----------------------------------------------------------------------------- 242 # 1. mu 243 xlim = range( c( mu ) ) 244 histInfo = plotPost( mu , cex.lab = 1.75 , 245 showCurve=showCurve , compVal=compVal , 246 xlab=bquote(mu) , main=paste("Mean") , 247 col="skyblue" ) 248 249 #----------------------------------------------------------------------------- 250 # 2. Plot data y and smattering of posterior predictive curves: 251 # Data histogram: 252 histInfo = hist( y , prob=TRUE , xlim=xLim , ylim=c(0,yMax) , breaks=xBreaks, 253 col="pink" , border="white" , xlab="y" , ylab="" , 254 yaxt="n" , cex.lab=1.5 , 255 main=paste("Data w. Post. Pred.") ) 256 # Posterior predictive probabilities of outcomes: 257 outProb=matrix(0,nrow=chainLength,ncol=max(y)) 258 for ( stepIdx in 1:chainLength ) { 259 threshCumProb = ( pnorm( ( mcmcMat[ stepIdx , 260 paste("thresh[",i,",",1:(max(y)-1),"]",sep="") ] 261 - mu[stepIdx] ) 262 / sigma[stepIdx] ) ) 263 outProb[stepIdx,] = c(threshCumProb,1) - c(0,threshCumProb) 264 } 265 outHdi = apply( outProb , 2 , HDIofMCMC ) 266 outMean = apply( outProb , 2 , median , na.rm=TRUE ) 267 show(outMean) 268 points( x=1:max(y) , y=outMean , pch=19 , cex=2 , col="skyblue" ) 269 segments( x0=1:max(y) , y0=outHdi[1,] , 270 x1=1:max(y) , y1=outHdi[2,] , lwd=4 , col="skyblue" ) 271 # annotate N: 272 text( max(xVec) , yMax , bquote(N==.(length(y))) , adj=c(1.1,1.1) ) 273 274 #----------------------------------------------------------------------------- 275 # 3. sigma 276 xlim=range( c( sigma ) ) 277 histInfo = plotPost( sigma , xlim=xlim , cex.lab = 1.75 , 278 showCurve=showCurve , 279 xlab=bquote(sigma) , 280 main=paste("Std. Dev.") , 281 col="skyblue" ) 282 283 #----------------------------------------------------------------------------- 284 # 4. effect size. 285 effectSize = ( mu - compVal ) / sigma 286 histInfo = plotPost( effectSize , compVal=0.0 , # ROPE=RopeEff , 287 showCurve=showCurve , 288 xlab=bquote( ( mu - .(compVal) ) / sigma ) , 289 cex.lab=1.75 , main="Effect Size" , 290 col="skyblue" ) 291 #----------------------------------------------------------------------------- 292 # 5. Thresholds: 293 threshCols = grep(paste0("^thresh\\[", i), colnames(mcmcMat), value=TRUE) 294 threshMean = rowMeans( mcmcMat[,threshCols] ) 295 xLim = range(mcmcMat[,threshCols]) 296 nPtToPlot = 500 297 plotIdx = floor(seq(1,nrow(mcmcMat),length=nPtToPlot)) 298 plot( mcmcMat[plotIdx,threshCols[1]] , threshMean[plotIdx] , col="skyblue" , 299 xlim=xLim , xlab="Threshold" , ylab="Mean Threshold" ) 300 abline(v=mean(mcmcMat[plotIdx,threshCols[1]]),lty="dashed",col="skyblue") 301 for ( jj in 2:length(threshCols) ) { 302 points( mcmcMat[plotIdx,threshCols[jj]] , threshMean[plotIdx] , col="skyblue" ) 303 abline(v=mean(mcmcMat[plotIdx,threshCols[jj]]),lty="dashed",col="skyblue") 304 } 305 306 #----------------------------------------------------------------------------- 307 # 6. (Is blank) 308 309 if ( !is.null(saveName) ) { 310 saveGraph( file=paste(saveName,"Q",i,"Post",sep=""), type=saveType) 311 } 312 } 313 } 314 315 #===============================================================================