### R code from vignette source 'lecture7_2015.Rnw' ### Encoding: UTF-8 ################################################### ### code chunk number 1: opts ################################################### options(useFancyQuotes=FALSE) ################################################### ### code chunk number 2: mcmcmod ################################################### library(foreign) dat <- read.dta("http://www.quantoid.net/files/essex/dr.dta") library(MCMCpack) inits1 <- runif(6,-2,2) inits2 <- runif(6,-2,2) chain1 <- MCMCregress(rep_mean ~ voice_mean*veto_mean + cwar + iwar, data=dat, burnin = 100000, thin=10, mcmc=100000, beta.start=rep(0,6), seed=1234) chain2 <- MCMCregress(rep_mean ~ voice_mean*veto_mean + cwar + iwar, data=dat, burnin = 100000, thin=10, mcmc=100000, beta.start=inits2, seed=5678) chains <- mcmc.list(chain1, chain2) ################################################### ### code chunk number 3: sums ################################################### summary(chains) ################################################### ### code chunk number 4: pc1 ################################################### plot(chains[,2], main="Voice") ################################################### ### code chunk number 5: acf1 ################################################### acf(chains[[1]][,1]) ################################################### ### code chunk number 6: gew ################################################### geweke.diag(chains[[1]]) ################################################### ### code chunk number 7: heidel ################################################### heidel.diag(chains[[1]]) ################################################### ### code chunk number 8: raft ################################################### raftery.diag(chains[[1]]) ################################################### ### code chunk number 9: gelman ################################################### gelman.diag(chains) ################################################### ### code chunk number 10: hpd1 ################################################### library(runjags) both <- combine.mcmc(chains) HPDinterval(both) ################################################### ### code chunk number 11: bayesFA (eval = FALSE) ################################################### ## dr <- read.dta("http://www.quantoid.net/files/essex/demrep.dta") ## cons <- list(xconst = c(2,0), xconst = list(1, "+"), polconiii = c(2,0), ## lgates = c(2,0), log_checks=c(2,0), disap = c(1,0), kill = c(1,0), ## tort=c(1,0), tort=list(2,"+")) ## X <- as.data.frame(scale(dr[which(dr$year == 2000),-(1:3)])) ## linits1 <- linits2 <- matrix(0, ncol=2, nrow=8) ## linits1[cbind(c(2,3,4,5,5,6,8), c(1,1,1,1,2,2,2) )] <- runif(7, 0,1) ## linits2[cbind(c(2,3,4,5,5,6,8), c(1,1,1,1,2,2,2) )] <- runif(7, 0,1) ## ## fa1 <- MCMCfactanal(~., X, factors=2, burnin=100000, lambda.start = linits1, ## mcmc=50000, thin=5, lambda.constraints=cons, seed=2002) ## fa2 <- MCMCfactanal(~., X, factors=2, burnin=100000, lambda.start = linits2, ## mcmc=50000, thin=5, lambda.constraints=cons, seed=1001) ## chains <- mcmc.list(fa1, fa2) ################################################### ### code chunk number 12: loadfa ################################################### load("fares.rda") ################################################### ### code chunk number 13: readdat ################################################### fh <- read.csv("http://quantoid.net/files/essex/jpr_replication.csv") fh06 <- fh[which(fh$year == 2006), ] X <- fh06[, c("A", "B", "C")] cons <- list(A=list(1, "+")) linits1 <- matrix(runif(3,0,1), ncol=1) linits2 <- matrix(runif(3,0,1), ncol=1) fh1a <- MCMCfactanal(~., data=as.data.frame(scale(X)), factors=1, burnin=25000, mcmc=50000, thin=5, lambda.start=linits1, lambda.constraint = cons, seed = 59302) fh2a <- MCMCfactanal(~., data=as.data.frame(scale(X)), factors=1, burnin=25000, mcmc=50000, thin=5, lambda.start=linits2, lambda.constraint = cons, seed = 10683) chains <- mcmc.list(fh1a, fh2a) fh2b <- MCMCfactanal(~., data=as.data.frame(scale(X)), factors=1, burnin=25000, mcmc=50000, thin=5, lambda.start=linits2, lambda.constraint = cons, store.scores=T, seed = 10683) fa1 <- fa(scale(X), nfac=1, fm="pa", scores=TRUE) cor(fa1$scores, colMeans(fh2b)[-(1:6)]) fh06 <- fh[which(fh$year == 2006), ] X <- fh06[, c("A", "B", "C", "D", "E", "F")] cons <- list(A=list(1, "+"), A=c(2,0), F=c(1,0), F=list(2, "+")) linits1 <- matrix(runif(12,0,1), ncol=2) linits2 <- matrix(runif(12,0,1), ncol=2) fh1a2 <- MCMCfactanal(~., data=as.data.frame(scale(X)), factors=2, burnin=25000, mcmc=50000, thin=5, lambda.start=linits1, lambda.constraint = cons, seed = 59302) fh2a2 <- MCMCfactanal(~., data=as.data.frame(scale(X)), factors=2, burnin=25000, mcmc=50000, thin=5, lambda.start=linits2, lambda.constraint = cons, seed = 10683) chains2 <- mcmc.list(fh1a2, fh2a2) fh06 <- fh[which(fh$year == 2006), ] X <- fh06[, c("A", "B", "C", "D", "E", "F")] cons <- list( A=c(1,1), A=c(2,0), F=c(1,0), F=c(2,1) ) linits1 <- matrix(runif(12,0,1), ncol=2) linits2 <- matrix(runif(12,0,1), ncol=2) linits1[1,1] <- NA linits1[1,2] <- NA linits1[6,1] <- NA linits1[6,2] <- NA linits2[1,1] <- NA linits2[1,2] <- NA linits2[6,1] <- NA linits2[6,2] <- NA fh1b2 <- MCMCfactanal(~., data=as.data.frame(scale(X)), factors=2, burnin=10000, mcmc=25000, thin=5, L0 = 1, a0=1, b0=1, lambda.start=linits1, lambda.constraint = cons, seed = 59302) fh2b2 <- MCMCfactanal(~., data=as.data.frame(scale(X)), factors=2, burnin=10000, mcmc=25000, thin=5, L0=1, a0=1, b0=1, lambda.start=linits2, lambda.constraint = cons, seed = 10683) chains2b <- mcmc.list(fh2b2, fh1b2) s <- seq(.001,.999, length=250) plot(s, densigamma(s, 1, 1)) fh2s <- MCMCfactanal(~., data=as.data.frame(scale(X)), factors=2, burnin=10000, mcmc=10000, thin=2, L0=1, a0=1, b0=1, lambda.start=linits2, lambda.constraint = cons, store.scores=TRUE, seed = 10683) fa2 <- fa(scale(X), nfac=2, fm="pa", scores=TRUE, rotate="promax") bfascores <- matrix(colMeans(fh2s)[grep("phi", colnames(fh2s))], ncol=2, byrow=T) cor(fa2$scores, bfascores) rownames(bfascores) <- fh06$Country bfascores <- as.data.frame(bfascores) bfascores <- bfascores[order(bfascores[,1]), ] bfascores$country <- factor(1:nrow(bfascores), labels=rownames(bfascores)) scores <- fh2s[,grep("phi", colnames(fh2s))] cis <- t(apply(scores, 2, quantile, c(.5,.025,.975))) scores <- cis[seq(1, 386, by=2), ] rownames(scores) <- fh06$Country scores <- scores[order(scores[,1]), ] scores <- as.data.frame(scores) scores$country <- factor(1:nrow(scores), labels=rownames(scores)) names(scores) <- c("med", "lower", "upper", "country") xl <- range(scores[,c("lower", "upper")])*1.1 pdf("dotplot1.pdf", height=8, width=5) dotplot(country ~ V1, data=bfascores[1:100, ], lty=0, scales=list(y=list(cex=.5)), col="black", cex=.5, xlim=range(bfascores$V1)*1.1) dev.off() pdf("dotplot2.pdf", height=8, width=5) dotplot(country ~ V1, data=bfascores[101:193, ], lty=0, scales=list(y=list(cex=.5)), col="black", cex=.5, xlim=range(bfascores$V1)*1.1) dev.off() pdf("dotplot1.pdf", height=8, width=5) dotplot(country ~ med, data=scores[1:100, ], lty=0, scales=list(y=list(cex=.5)), col="black", cex=.5, xlim=xl) trellis.focus() panel.segments(scores$lower[1:100], 1:100, scores$upper[1:100], 1:100) trellis.unfocus() dev.off() pdf("dotplot2.pdf", height=8, width=5) dotplot(country ~ V1, data=bfascores[101:193, ], lty=0, scales=list(y=list(cex=.5)), col="black", cex=.5, xlim=range(bfascores$V1)*1.1) dev.off()