source("04jspid.txt") scps.data <- list("N", "M", "R", "county", "Y", "imal", "iblk", "icvs", "icvm", "ictr", "icdr", "icpr", "icjs", "ipri", "idet", "irev", "itri", "cidx", "cunr", "cblp", "ccon", "csth", "cgui") scps.inits <- function(){ list(theta=rnorm(26), eta.raw=rnorm(86), b=structure(.Data=rnorm(507), .Dim=c(39,13)), gamma=structure(.Data=diag(rlnorm(13)), .Dim=c(13,13)), xi=rnorm(13), imal=replace(replace(imal, which(is.na(imal)), 0), which(!is.na(imal)), NA), iblk=replace(replace(iblk, which(is.na(iblk)), 0), which(!is.na(iblk)), NA), icjs=replace(replace(icjs, which(is.na(icjs)), 0), which(!is.na(icjs)), NA), ipri=replace(replace(ipri, which(is.na(ipri)), 0), which(!is.na(ipri)), NA), idet=replace(replace(idet, which(is.na(idet)), 0), which(!is.na(idet)), NA)) } scps.parameters <- c("beta", "e.y", "e.b", "theta", "eta") scps.sims <- bugs(scps.data, scps.inits, scps.parameters, "04jspim.txt", n.chains=3, n.iter=10000, n.thin=10) theta[1] 1.7 0.0 1.6 1.7 1.7 1.7 1.8 1.0 1500 theta[2] 0.6 0.2 0.2 0.5 0.6 0.7 1.0 1.0 950 theta[3] -0.3 0.1 -0.5 -0.4 -0.3 -0.3 -0.2 1.0 770 theta[4] -0.5 0.1 -0.6 -0.5 -0.5 -0.4 -0.3 1.0 1500 theta[5] 0.6 0.1 0.4 0.5 0.6 0.7 0.9 1.0 1500 theta[6] 0.0 0.0 0.0 0.0 0.0 0.1 0.1 1.0 930 theta[7] 0.4 0.1 0.2 0.3 0.4 0.5 0.7 1.0 1500 theta[8] 0.6 0.1 0.5 0.5 0.6 0.6 0.7 1.0 1500 theta[9] -0.1 0.1 -0.2 -0.1 -0.1 -0.1 0.0 1.0 1500 theta[10] 0.1 0.1 0.0 0.1 0.1 0.2 0.3 1.0 1500 theta[11] 0.4 0.1 0.2 0.3 0.4 0.4 0.6 1.0 1500 theta[12] -0.4 0.0 -0.5 -0.4 -0.4 -0.4 -0.3 1.0 1500 theta[13] -0.2 0.1 -0.4 -0.3 -0.2 -0.2 -0.1 1.0 1500 theta[14] 0.2 0.1 0.0 0.1 0.2 0.2 0.3 1.0 1500 theta[15] 0.1 0.1 0.0 0.1 0.1 0.1 0.2 1.0 1500 theta[16] -0.4 0.1 -0.6 -0.4 -0.4 -0.3 -0.2 1.0 1500 theta[17] -0.8 0.0 -0.8 -0.8 -0.8 -0.8 -0.7 1.0 710 theta[18] -0.5 0.1 -0.7 -0.6 -0.5 -0.4 -0.3 1.0 1500 theta[19] -0.2 0.1 -0.4 -0.3 -0.2 -0.2 -0.1 1.0 740 theta[20] -0.5 0.0 -0.6 -0.6 -0.5 -0.5 -0.5 1.0 1300 theta[21] 1.7 0.1 1.5 1.6 1.7 1.8 2.0 1.0 1500 theta[22] 0.8 0.1 0.7 0.8 0.8 0.9 1.0 1.0 1500 theta[23] 0.4 0.1 0.3 0.4 0.4 0.5 0.6 1.0 1100 theta[24] 0.2 0.1 0.1 0.2 0.2 0.3 0.3 1.0 1100 theta[25] 0.5 0.1 0.3 0.4 0.5 0.5 0.6 1.0 1500 theta[26] 0.2 0.1 0.0 0.2 0.2 0.3 0.4 1.0 1500 eta[1] -5.2 0.4 -5.9 -5.4 -5.1 -4.9 -4.5 1.0 330 eta[2] 0.4 0.3 -0.3 0.1 0.3 0.6 1.0 1.0 51 eta[3] -0.6 0.4 -1.3 -0.9 -0.6 -0.4 0.1 1.0 200 eta[4] 0.0 0.3 -0.6 -0.2 0.0 0.3 0.7 1.0 120 eta[5] 0.6 0.3 -0.1 0.3 0.6 0.8 1.2 1.0 47 eta[6] -0.7 0.6 -1.9 -1.0 -0.7 -0.3 0.5 1.0 550 eta[7] 0.1 0.6 -0.9 -0.2 0.1 0.5 1.3 1.0 240 eta[8] 0.5 0.1 0.3 0.4 0.5 0.6 0.7 1.0 89 eta[9] 0.0 0.1 -0.1 0.0 0.0 0.0 0.2 1.0 470 eta[10] 0.0 0.1 -0.2 -0.1 0.0 0.0 0.2 1.1 27 eta[11] 0.0 0.1 -0.1 0.0 0.0 0.1 0.2 1.1 23 eta[12] 0.0 0.1 -0.1 0.0 0.0 0.1 0.2 1.0 340 eta[13] 0.1 0.1 -0.1 0.0 0.1 0.1 0.3 1.1 41 eta[14] 0.0 0.1 -0.2 -0.1 0.0 0.1 0.2 1.0 410 eta[15] 0.2 0.1 0.0 0.2 0.2 0.3 0.5 1.0 60 eta[16] 0.0 0.1 -0.2 -0.1 0.0 0.1 0.2 1.0 76 eta[17] 0.1 0.1 -0.1 0.0 0.1 0.1 0.3 1.0 72 eta[18] 0.0 0.1 -0.2 -0.1 0.0 0.1 0.2 1.0 140 eta[19] 0.0 0.1 -0.3 -0.1 0.0 0.1 0.3 1.0 620 eta[20] 2.6 0.3 2.1 2.4 2.6 2.7 3.1 1.0 47 eta[21] -0.1 0.2 -0.6 -0.3 -0.1 0.0 0.3 1.0 140 eta[22] -0.1 0.3 -0.6 -0.3 -0.1 0.1 0.6 1.0 570 eta[23] 0.5 0.2 0.0 0.3 0.5 0.6 1.0 1.0 420 eta[24] 0.2 0.3 -0.3 0.0 0.2 0.4 0.7 1.1 34 eta[25] 0.0 0.3 -0.8 -0.2 0.0 0.2 0.6 1.0 130 eta[26] 0.0 0.4 -0.8 -0.2 0.1 0.3 0.7 1.1 38 eta[27] 1.6 0.2 1.2 1.4 1.6 1.7 2.0 1.0 130 eta[28] -0.3 0.2 -0.7 -0.5 -0.3 -0.2 0.0 1.0 1500 eta[29] 0.2 0.2 -0.2 0.1 0.2 0.3 0.6 1.0 260 eta[30] 0.4 0.2 0.0 0.3 0.4 0.5 0.8 1.0 1500 eta[31] 0.1 0.2 -0.3 0.0 0.1 0.2 0.4 1.1 32 eta[32] 0.5 0.3 -0.1 0.3 0.5 0.7 1.2 1.0 480 eta[33] -0.1 0.3 -0.7 -0.3 -0.1 0.2 0.6 1.0 1500 eta[34] 1.5 0.2 1.1 1.4 1.5 1.7 2.0 1.1 37 eta[35] 0.0 0.2 -0.4 -0.2 0.0 0.1 0.4 1.1 41 eta[36] -0.2 0.2 -0.7 -0.4 -0.2 -0.1 0.3 1.0 140 eta[37] 0.0 0.2 -0.4 -0.1 0.0 0.2 0.4 1.0 110 eta[38] 0.0 0.2 -0.4 -0.2 0.0 0.1 0.5 1.0 48 eta[39] -0.1 0.3 -0.8 -0.3 -0.1 0.1 0.5 1.0 230 eta[40] 0.5 0.4 -0.2 0.2 0.5 0.7 1.2 1.0 79 eta[41] 0.4 0.2 0.0 0.3 0.4 0.6 0.9 1.1 24 eta[42] -0.3 0.2 -0.8 -0.5 -0.3 -0.2 0.1 1.0 100 eta[43] 0.2 0.3 -0.3 0.0 0.2 0.4 0.7 1.1 27 eta[44] 0.4 0.2 0.0 0.2 0.4 0.5 0.9 1.0 130 eta[45] 0.3 0.2 -0.1 0.2 0.3 0.5 0.8 1.1 42 eta[46] 0.1 0.3 -0.5 -0.1 0.1 0.3 0.7 1.1 45 eta[47] 0.2 0.4 -0.5 0.0 0.2 0.5 1.0 1.0 97 eta[48] 1.4 0.2 1.0 1.3 1.4 1.6 1.9 1.0 66 eta[49] -0.2 0.2 -0.6 -0.4 -0.2 -0.1 0.1 1.0 51 eta[50] 0.0 0.2 -0.4 -0.1 0.0 0.1 0.4 1.1 24 eta[51] 0.3 0.2 -0.1 0.1 0.3 0.4 0.6 1.0 110 eta[52] -0.1 0.2 -0.5 -0.3 -0.1 0.0 0.3 1.1 29 eta[53] -0.3 0.3 -0.9 -0.5 -0.3 -0.1 0.3 1.0 140 eta[54] -0.7 0.3 -1.4 -0.9 -0.7 -0.5 0.0 1.0 1500 eta[55] 0.8 0.1 0.6 0.7 0.8 0.9 1.0 1.1 41 eta[56] 0.0 0.1 -0.1 0.0 0.0 0.1 0.2 1.0 160 eta[57] 0.0 0.1 -0.2 -0.1 0.0 0.1 0.2 1.0 390 eta[58] -0.1 0.1 -0.3 -0.1 -0.1 0.0 0.1 1.1 54 eta[59] 0.1 0.1 -0.1 0.0 0.1 0.1 0.3 1.1 35 eta[60] 0.1 0.2 -0.2 0.0 0.1 0.2 0.4 1.1 250 eta[61] 1.7 0.1 1.5 1.7 1.7 1.8 2.0 1.0 100 eta[62] -0.1 0.1 -0.3 -0.2 -0.1 0.0 0.2 1.1 41 eta[63] 0.3 0.1 0.0 0.2 0.3 0.4 0.6 1.0 150 eta[64] -0.3 0.1 -0.5 -0.4 -0.3 -0.2 -0.1 1.0 84 eta[65] -0.1 0.1 -0.3 -0.1 -0.1 0.0 0.2 1.1 44 eta[66] -0.2 0.2 -0.6 -0.3 -0.2 0.0 0.2 1.0 870 eta[67] -0.1 0.2 -0.5 -0.3 -0.1 0.0 0.3 1.0 62 eta[68] 1.9 0.1 1.7 1.8 1.9 2.0 2.2 1.0 970 eta[69] 0.1 0.2 -0.2 0.0 0.1 0.2 0.4 1.1 18 eta[70] -0.2 0.2 -0.5 -0.3 -0.1 0.0 0.2 1.0 130 eta[71] -0.1 0.1 -0.3 -0.2 -0.1 0.0 0.2 1.0 1500 eta[72] 0.1 0.2 -0.3 0.0 0.1 0.2 0.4 1.0 78 eta[73] -0.5 0.2 -1.0 -0.7 -0.5 -0.4 -0.1 1.0 120 eta[74] 1.3 0.2 1.0 1.2 1.3 1.4 1.7 1.0 1500 eta[75] 0.1 0.2 -0.2 0.0 0.1 0.2 0.5 1.2 18 eta[76] 0.2 0.2 -0.1 0.1 0.2 0.3 0.7 1.0 62 eta[77] -0.1 0.2 -0.4 -0.2 -0.1 0.0 0.3 1.0 1500 eta[78] 0.2 0.2 -0.2 0.1 0.2 0.4 0.7 1.1 32 eta[79] -0.2 0.3 -0.7 -0.3 -0.2 0.0 0.3 1.1 64 eta[80] 0.7 0.2 0.2 0.5 0.7 0.8 1.1 1.0 120 eta[81] -0.2 0.2 -0.6 -0.3 -0.2 -0.1 0.1 1.0 61 eta[82] 0.0 0.2 -0.5 -0.2 0.0 0.1 0.4 1.0 540 eta[83] -0.2 0.2 -0.6 -0.3 -0.2 0.0 0.2 1.0 120 eta[84] 0.1 0.2 -0.3 0.0 0.1 0.2 0.4 1.1 37 eta[85] 0.2 0.3 -0.3 0.0 0.2 0.3 0.8 1.0 88 eta[86] 0.1 0.3 -0.4 -0.1 0.0 0.2 0.6 1.0 1200 deviance 51990.8 41.7 51910.0 51960.0 51990.0 52020.0 52070.0 1.0 640 pD = 866.4 and DIC = 52857.2 (using the rule, pD = var(deviance)/2) For each parameter, n.eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor (at convergence, Rhat=1). DIC is an estimate of expected predictive error (lower deviance is better). rsquared.y <- 1 - mean (apply (e.y, 1, var)) / var (Y) lambda.y <- 1 - var (apply (e.y, 2, mean)) / mean (apply (e.y, 1, var)) rsquared.b <- lambda.b <- rep(NA,13) for (k in 1:13) { rsquared.b[k] <- 1 - mean (apply (e.b[,,k], 1, var)) / mean (apply (beta[,,k], 1, var)) lambda.b[k] <- 1 - var (apply (e.b[,,k], 2, mean)) / mean (apply (e.b[,,k], 1, var)) } print(round(c(rsquared.y,rsquared.b), 2)) #0.40 0.38 0.61 0.82 0.43 0.70 0.19 0.44 0.38 0.43 0.77 0.46 0.35 0.58 print(round(c(lambda.y,lambda.b), 2)) #0.06 0.34 0.92 0.93 0.74 0.86 0.56 0.62 0.58 0.80 0.86 0.68 0.67 0.87 prismeans <- apply(eta, 2, mean) thetameans <- apply(theta, 2, mean) alpha <- e.b alphameans <- apply(alpha, c(2, 3), mean) sample.indices <- sample(dim(eta)[1], size=100) eta <- eta[sample.indices,] alpha <- alpha[sample.indices,,] dump(c("eta", "alpha", "prismeans", "alphameans", "thetameans"), file="resultsh.txt")