### Code for the correlation plot in Zheng et al (2006) ### The "How many prisoners do you know" paper. ### [IMPORTANT] You may need to re-size your graphics window ### after the plot is done to get better effects. ############### ### Step 0: get plot data attach(sims.all, name='sims.obj') n.keep<-n.iter/n.thin result.alpha<-array(NA, c(m.iter, data.n, n.keep/2)) result.beta<-array(NA, c(m.iter, data.j, n.keep/2)) result.omega<-array(NA, c(m.iter, data.j, n.keep/2)) for (i in 1:data.n) result.alpha[,i,] <- t(sims[-(1:(n.keep/2)),,i]) for (i in 1:data.j) {result.beta[,i,] <- t(sims[-(1:(n.keep/2)),,data.n+i]) result.omega[,i,] <- t(sims[-(1:(n.keep/2)),,data.n+data.j+i])} y.pred<-exp(outer(result.alpha[1,,105], result.beta[1,,105], '+')) y.diff<-sqrt(y)-sqrt(y.pred) y.std<-y.diff for(j in 1:data.j){ y.std[,j]<-(y.diff[,j]-mean(y.diff[,j], na.rm=T))/sqrt(var(y.diff[,j], use='complete')) } ### Step 1: prepare correlation matrix z.plot<-cor(sqrt(y)-sqrt(y.pred), use='pairwise.complete') ### Step 2: calculate order of the rows/columns ### it can be specified manually or computed using hierarchical clustering. z.order<-c(11,9,7,5,3,1,12,10,8,6,4,2) ## first names #z.other<-z.other[hclust(dist(t(y.std[,z.other])))\$order] z.other<-c(32,31,28,27,30,29,25,26,13,20,24,18,16,21,17,19,23,22,14,15) ## other category z.order<-c(z.other, z.order) ### Step 3: order the correlation matrix according to z.order z.plot<-z.plot[z.order, z.order] ### Step 4: mask out the upper half for(i in 1:32) for(j in i:32) z.plot[i,j]<-0 ## in this project, we do not pay attention to negative correlation ### Step 5: get the range of the correlation values max.cor<-0.7 min.cor<-min(z.plot) ### Step 6: graphics layout layout(matrix(c(1,2), 1, 2, byrow=TRUE), c(10,1)) par(mar=c(0.1, 0.1, 2, 0.1), pty='s') ### break values for image function ### here the number of levels is set to 8 z.breaks<-c(min(z.plot), seq(0, 0.7, 0.1)) z.colors<-gray(8:1/8) image(1:32, 1:32, z.plot, xaxt='n', yaxt='n', col=z.colors, breaks=z.breaks, xlim=c(-6, 33.5), ylim=c(-6, 33.5), xaxt='n', yaxt='n', bty='n', xlab='', ylab='') ### Step 7: plot legends ### category names text(1:32, 1:32, category.short[z.order], adj=1, cex=0.9) ### boxes rect(26.5, 26.5, 32.5, 32.5, lty=2) text(26, 32.5, 'female \n names', cex=1, adj=1, font=2) rect(20.5, 20.5, 26.5, 26.5, lty=2) text(20, 26.5, 'male \n names', cex=1, adj=1, font=2) rect(8.5, 8.5, 0.5, 0.5, lty=2) text(0, 8.5, 'negative \n experiences', adj=1, cex=1, font=2) ### the color legend par(mar=c(4, 0.1, 2 ,0.1), pty='m') plot(c(0.2,1), c(min(z.breaks), max(z.breaks)), type='n', bty='n', xlab='', ylab='', xaxt='n', yaxt='n') options(digits=1) for(i in 2:(length(z.breaks))){ rect(0.5, z.breaks[i-1], 1, z.breaks[i], col=z.colors[i-1]) text(0.45, z.breaks[i-1], format(round(z.breaks[i-1],1)), cex=1, adj=1) } rect(0.5, z.breaks[length(z.breaks)], 1, z.breaks[length(z.breaks)], col=z.colors[length(z.colors)]) detach('sims.obj') #######################################