##################################
# Splus code for 2004 CSDA paper #
##################################

################
# Read in data #
################

prison98 <- read.table("prison98.txt", header = T, row.names = NULL)
attach(prison98)

IMAL.na <- is.na(IMAL)
IBLK.na <- is.na(IBLK)
ICJS.na <- is.na(ICJS)
IPRI.na <- is.na(IPRI)
IDET.na <- is.na(IDET)

meanidx <- 587.333333
sdidx <- 219.840335
meanunr <- 4.428205
sdunr <- 1.759478
meanblp <- 18.879487
sdblp <- 12.380719
meancon <- 38.225641
sdcon <- 13.256373

Cidx <- meanidx + sdidx * CIDX
Cunr <- meanunr + sdunr * CUNR
Cblp <- meanblp + sdblp * CBLP
Ccon <- meancon + sdcon * CCON

idx <- unr <- blp <- con <- NULL
for (i in 1:39) {
	idx <- c(idx, CIDX[COUNTY==i][1])
	unr <- c(unr, CUNR[COUNTY==i][1])
	blp <- c(blp, CBLP[COUNTY==i][1])
	con <- c(con, CCON[COUNTY==i][1])
}
idx <- meanidx + sdidx * idx
unr <- meanunr + sdunr * unr
blp <- meanblp + sdblp * blp
con <- meancon + sdcon * con

#############################################################
# Hierarchical results, posterior means and 100 simulations #
#############################################################

source("resultsh.txt")
eta.samp <- t(eta)
alpha.samp <- t(array(aperm(alpha,c(1,3,2)),c(100,507)))

round(apply(alphameans, 2, mean), 2) # means over counties
# -0.02 0.00 0.00 0.02 0.01 0.02 0.02 0.03 0.01 0.00 0.00 -0.01 0.01

round(apply(array(apply(alpha, 1, var),c(13,13,100)),c(1,2),mean),2)
# var-cov matrices over counties
# 1.52-0.02-0.01 0.00-0.01 0.19 0.10-0.13-0.01-0.01-0.12-0.02 0.09
#-0.02 0.02 0.00-0.01-0.01-0.04-0.04-0.04 0.00 0.00 0.01 0.02-0.02
#-0.01 0.00 0.03-0.03 0.00-0.01 0.01 0.01 0.00 0.01 0.01 0.00 0.00
# 0.00-0.01-0.03 0.47 0.02 0.29 0.10 0.11-0.04-0.01-0.16-0.14-0.02
#-0.01-0.01 0.00 0.02 0.19 0.15 0.06 0.13 0.00-0.02-0.04-0.06 0.03
# 0.19-0.04-0.01 0.29 0.15 0.71 0.38 0.41 0.01-0.03-0.23-0.28 0.10
# 0.10-0.04 0.01 0.10 0.06 0.38 0.53 0.39 0.05 0.02-0.08-0.21 0.13
#-0.13-0.04 0.01 0.11 0.13 0.41 0.39 0.56 0.06 0.02-0.10-0.24 0.11
#-0.01 0.00 0.00-0.04 0.00 0.01 0.05 0.06 0.05 0.01 0.01-0.02 0.02
#-0.01 0.00 0.01-0.01-0.02-0.03 0.02 0.02 0.01 0.07 0.01-0.01-0.01
#-0.12 0.01 0.01-0.16-0.04-0.23-0.08-0.10 0.01 0.01 0.19 0.15 0.00
#-0.02 0.02 0.00-0.14-0.06-0.28-0.21-0.24-0.02-0.01 0.15 0.24-0.04
# 0.09-0.02 0.00-0.02 0.03 0.10 0.13 0.11 0.02-0.01 0.00-0.04 0.21

alphamean2 <- apply(alpha, c(2, 3), mean) # means over samples

trellis.device()

for (i in 1:13) {
	qqplot <- qqmath(~alphamean2[,i],
   	    prepanel = prepanel.qqmathline,
   	    panel = function(x, y) {
   	      panel.grid()
   	      panel.qqmathline(y, distribution = qnorm)
   	      panel.qqmath(x, y)
   	    },
	)
	print(qqplot)
}

for(i in 1:13) {
	plot(idx, alphamean2[, i])
}

for(i in 1:13) {
	plot(unr, alphamean2[, i])
}

for(i in 1:13) {
	plot(blp, alphamean2[, i])
}

for(i in 1:13) {
	plot(con, alphamean2[, i])
}

###################################################
# Missing data imputation from hierarchical model #
###################################################

IMAL.im <- 1/(1+exp(-thetameans[1] - 
	thetameans[2]*ICVS[IMAL.na] -
	thetameans[3]*ICDR[IMAL.na] -
	thetameans[4]*ICPR[IMAL.na] -
	thetameans[5]*ITRI[IMAL.na]))
IMAL.h <- replace(IMAL, names(IMAL[IMAL.na]), IMAL.im)

IBLK.im <- 1/(1+exp(-thetameans[6] - 
	thetameans[7]*ICVS[IBLK.na] -
	thetameans[8]*ICTR[IBLK.na] -
	thetameans[9]*ICPR[IBLK.na] -
	thetameans[10]*IREV[IBLK.na] -
	thetameans[11]*ITRI[IBLK.na]))
IBLK.h <- replace(IBLK, names(IBLK[IBLK.na]), IBLK.im)

ICJS.im <- 1/(1+exp(-thetameans[12] - 
	thetameans[13]*ICVM[ICJS.na] -
	thetameans[14]*ICDR[ICJS.na] -
	thetameans[15]*ICPR[ICJS.na] -
	thetameans[16]*ITRI[ICJS.na]))
ICJS.h <- replace(ICJS, names(ICJS[ICJS.na]), ICJS.im)

IPRI.im <- 1/(1+exp(-thetameans[17] - 
	thetameans[18]*ICVM[IPRI.na] -
	thetameans[19]*IREV[IPRI.na]))
IPRI.h <- replace(IPRI, names(IPRI[IPRI.na]), IPRI.im)

IDET.im <- 1/(1+exp(-thetameans[20] - 
	thetameans[21]*ICVS[IDET.na] -
	thetameans[22]*ICVM[IDET.na] -
	thetameans[23]*ICTR[IDET.na] -
	thetameans[24]*ICDR[IDET.na] -
	thetameans[25]*ICPR[IDET.na] -
	thetameans[26]*ITRI[IDET.na]))
IDET.h <- replace(IDET, names(IDET[IDET.na]), IDET.im)

##############################################
# Non-hierarchical results, 4000 simulations #
##############################################

boapris.n <- 
  read.table("resultsn.txt", header = T, row.names = 1, sep = ",")
attach(boapris.n)

beta <- NULL
for (i in 1:86)
  beta <- cbind(beta,
                eval(parse(text=paste("beta.", i, ".", sep=""))))
beta <- t(beta)

#######################################################
# Missing data imputation from non-hierarchical model #
#######################################################

IMAL.im <- 1/(1+exp(-mean(theta.1.) - 
	mean(theta.2.)*ICVS[IMAL.na] -
	mean(theta.3.)*ICDR[IMAL.na] -
	mean(theta.4.)*ICPR[IMAL.na] -
	mean(theta.5.)*ITRI[IMAL.na]))
IMAL.n <- replace(IMAL, names(IMAL[IMAL.na]), IMAL.im)

IBLK.im <- 1/(1+exp(-mean(theta.6.) - 
	mean(theta.7.)*ICVS[IBLK.na] -
	mean(theta.8.)*ICTR[IBLK.na] -
	mean(theta.9.)*ICPR[IBLK.na] -
	mean(theta.10.)*IREV[IBLK.na] -
	mean(theta.11.)*ITRI[IBLK.na]))
IBLK.n <- replace(IBLK, names(IBLK[IBLK.na]), IBLK.im)

ICJS.im <- 1/(1+exp(-mean(theta.12.) - 
	mean(theta.13.)*ICVM[ICJS.na] -
	mean(theta.14.)*ICDR[ICJS.na] -
	mean(theta.15.)*ICPR[ICJS.na] -
	mean(theta.16.)*ITRI[ICJS.na]))
ICJS.n <- replace(ICJS, names(ICJS[ICJS.na]), ICJS.im)

IPRI.im <- 1/(1+exp(-mean(theta.17.) - 
	mean(theta.18.)*ICVM[IPRI.na] -
	mean(theta.19.)*IREV[IPRI.na]))
IPRI.n <- replace(IPRI, names(IPRI[IPRI.na]), IPRI.im)

IDET.im <- 1/(1+exp(-mean(theta.20.) - 
	mean(theta.21.)*ICVS[IDET.na] -
	mean(theta.22.)*ICVM[IDET.na] -
	mean(theta.23.)*ICTR[IDET.na] -
	mean(theta.24.)*ICDR[IDET.na] -
	mean(theta.25.)*ICPR[IDET.na] -
	mean(theta.26.)*ITRI[IDET.na]))
IDET.n <- replace(IDET, names(IDET[IDET.na]), IDET.im)

####################
# Hierarchical fit #
####################

W.h <-
  cbind(1, CIDX, CUNR, CBLP, CCON, CSTH, CGUI, 
	IMAL.h, I(CIDX*IMAL.h), I(CUNR*IMAL.h), I(CBLP*IMAL.h),
        I(CCON*IMAL.h), I(CSTH*IMAL.h),  
	IBLK.h, I(CIDX*IBLK.h), I(CUNR*IBLK.h), I(CBLP*IBLK.h),
        I(CCON*IBLK.h), I(CSTH*IBLK.h),  
	ICVS, I(CIDX*ICVS), I(CUNR*ICVS), I(CBLP*ICVS), I(CCON*ICVS),
        I(CSTH*ICVS), I(CGUI*ICVS),  
	ICVM, I(CIDX*ICVM), I(CUNR*ICVM), I(CBLP*ICVM), I(CCON*ICVM),
        I(CSTH*ICVM), I(CGUI*ICVM),
        ICTR, I(CIDX*ICTR), I(CUNR*ICTR), I(CBLP*ICTR), I(CCON*ICTR),
        I(CSTH*ICTR), I(CGUI*ICTR),  
	ICDR, I(CIDX*ICDR), I(CUNR*ICDR), I(CBLP*ICDR), I(CCON*ICDR),
        I(CSTH*ICDR), I(CGUI*ICDR),  
	ICPR, I(CIDX*ICPR), I(CUNR*ICPR), I(CBLP*ICPR), I(CCON*ICPR),
        I(CSTH*ICPR), I(CGUI*ICPR),  
	ICJS.h, I(CIDX*ICJS.h), I(CUNR*ICJS.h), I(CBLP*ICJS.h),
        I(CCON*ICJS.h), I(CSTH*ICJS.h),  
	IPRI.h, I(CIDX*IPRI.h), I(CUNR*IPRI.h), I(CBLP*IPRI.h),
        I(CCON*IPRI.h), I(CSTH*IPRI.h), I(CGUI*IPRI.h),  
	IDET.h, I(CIDX*IDET.h), I(CUNR*IDET.h), I(CBLP*IDET.h),
        I(CCON*IDET.h), I(CSTH*IDET.h),  
	IREV, I(CIDX*IREV), I(CUNR*IREV), I(CBLP*IREV), I(CCON*IREV),
        I(CSTH*IREV),  
	ITRI, I(CIDX*ITRI), I(CUNR*ITRI), I(CBLP*ITRI), I(CCON*ITRI),
        I(CSTH*ITRI), I(CGUI*ITRI))

X1 <- cbind(1, IMAL.h, IBLK.h, ICVS, ICVM, ICTR, ICDR, ICPR, ICJS.h,
            IPRI.h, IDET.h, IREV, ITRI) 
X.h <- matrix(NA, nrow=8446, ncol=507)
for(i in 1:507) {
  X.h[,i] <- ifelse(COUNTY==(i-1)%/%13+1, X1[,(i-1)%%13+1], 0)
}

fit1.h <-
  glm(Y ~ CIDX+CUNR+CBLP+CCON+CSTH+CGUI+
      IMAL.h+I(CIDX*IMAL.h)+I(CUNR*IMAL.h)+I(CBLP*IMAL.h)+
      I(CCON*IMAL.h)+I(CSTH*IMAL.h)+
      IBLK.h+I(CIDX*IBLK.h)+I(CUNR*IBLK.h)+I(CBLP*IBLK.h)+
      I(CCON*IBLK.h)+I(CSTH*IBLK.h)+
      ICVS+I(CIDX*ICVS)+I(CUNR*ICVS)+I(CBLP*ICVS)+I(CCON*ICVS)+
      I(CSTH*ICVS)+I(CGUI*ICVS)+
      ICVM+I(CIDX*ICVM)+I(CUNR*ICVM)+I(CBLP*ICVM)+I(CCON*ICVM)+
      I(CSTH*ICVM)+I(CGUI*ICVM)+
      ICTR+I(CIDX*ICTR)+I(CUNR*ICTR)+I(CBLP*ICTR)+I(CCON*ICTR)+
      I(CSTH*ICTR)+I(CGUI*ICTR)+
      ICDR+I(CIDX*ICDR)+I(CUNR*ICDR)+I(CBLP*ICDR)+I(CCON*ICDR)+
      I(CSTH*ICDR)+I(CGUI*ICDR)+
      ICPR+I(CIDX*ICPR)+I(CUNR*ICPR)+I(CBLP*ICPR)+I(CCON*ICPR)+
      I(CSTH*ICPR)+I(CGUI*ICPR)+
      ICJS.h+I(CIDX*ICJS.h)+I(CUNR*ICJS.h)+I(CBLP*ICJS.h)+
      I(CCON*ICJS.h)+I(CSTH*ICJS.h)+
      IPRI.h+I(CIDX*IPRI.h)+I(CUNR*IPRI.h)+I(CBLP*IPRI.h)+
      I(CCON*IPRI.h)+I(CSTH*IPRI.h)+I(CGUI*IPRI.h)+ 
      IDET.h+I(CIDX*IDET.h)+I(CUNR*IDET.h)+I(CBLP*IDET.h)+
      I(CCON*IDET.h)+I(CSTH*IDET.h)+
      IREV+I(CIDX*IREV)+I(CUNR*IREV)+I(CBLP*IREV)+I(CCON*IREV)+
      I(CSTH*IREV)+
      ITRI+I(CIDX*ITRI)+I(CUNR*ITRI)+I(CBLP*ITRI)+I(CCON*ITRI)+
      I(CSTH*ITRI)+I(CGUI*ITRI),
      family = binomial)

fit1.h$linear.predictors <-
  (W.h %*% prismeans + X.h %*% as.vector(t(alphameans)))[,1]
fit1.h$fitted.values <- 1/(1 + exp(-fit1.h$linear.predictors))

linpred1.h <- W.h %*% eta.samp + X.h %*% alpha.samp
fitvals1.h <- list(yhat = 1/(1 + exp(-linpred1.h)))

########################
# Non-hierarchical fit #
########################

W.n <-
  cbind(1, CIDX, CUNR, CBLP, CCON, CSTH, CGUI, 
	IMAL.n, I(CIDX*IMAL.n), I(CUNR*IMAL.n), I(CBLP*IMAL.n),
        I(CCON*IMAL.n), I(CSTH*IMAL.n),  
	IBLK.n, I(CIDX*IBLK.n), I(CUNR*IBLK.n), I(CBLP*IBLK.n),
        I(CCON*IBLK.n), I(CSTH*IBLK.n),  
	ICVS, I(CIDX*ICVS), I(CUNR*ICVS), I(CBLP*ICVS), I(CCON*ICVS),
        I(CSTH*ICVS), I(CGUI*ICVS),  
	ICVM, I(CIDX*ICVM), I(CUNR*ICVM), I(CBLP*ICVM), I(CCON*ICVM),
        I(CSTH*ICVM), I(CGUI*ICVM),  
	ICTR, I(CIDX*ICTR), I(CUNR*ICTR), I(CBLP*ICTR), I(CCON*ICTR),
        I(CSTH*ICTR), I(CGUI*ICTR),  
	ICDR, I(CIDX*ICDR), I(CUNR*ICDR), I(CBLP*ICDR), I(CCON*ICDR),
        I(CSTH*ICDR), I(CGUI*ICDR),  
	ICPR, I(CIDX*ICPR), I(CUNR*ICPR), I(CBLP*ICPR), I(CCON*ICPR),
        I(CSTH*ICPR), I(CGUI*ICPR),  
	ICJS.n, I(CIDX*ICJS.n), I(CUNR*ICJS.n), I(CBLP*ICJS.n),
        I(CCON*ICJS.n), I(CSTH*ICJS.n),  
	IPRI.n, I(CIDX*IPRI.n), I(CUNR*IPRI.n), I(CBLP*IPRI.n),
        I(CCON*IPRI.n), I(CSTH*IPRI.n), I(CGUI*IPRI.n),  
	IDET.n, I(CIDX*IDET.n), I(CUNR*IDET.n), I(CBLP*IDET.n),
        I(CCON*IDET.n), I(CSTH*IDET.n),  
	IREV, I(CIDX*IREV), I(CUNR*IREV), I(CBLP*IREV), I(CCON*IREV),
        I(CSTH*IREV),  
	ITRI, I(CIDX*ITRI), I(CUNR*ITRI), I(CBLP*ITRI), I(CCON*ITRI),
        I(CSTH*ITRI), I(CGUI*ITRI))

fit2.n <-
  glm(Y ~ CIDX + CUNR + CBLP + CCON + CSTH + CGUI + 
      IMAL.n + I(CIDX*IMAL.n) + I(CUNR*IMAL.n) + I(CBLP*IMAL.n) + 
      I(CCON*IMAL.n) + I(CSTH*IMAL.n) + 
      IBLK.n + I(CIDX*IBLK.n) + I(CUNR*IBLK.n) + I(CBLP*IBLK.n) + 
      I(CCON*IBLK.n) + I(CSTH*IBLK.n) + 
      ICVS + I(CIDX*ICVS) + I(CUNR*ICVS) + I(CBLP*ICVS) + I(CCON*ICVS)
      +  I(CSTH*ICVS) + I(CGUI*ICVS) + 
      ICVM + I(CIDX*ICVM) + I(CUNR*ICVM) + I(CBLP*ICVM) + I(CCON*ICVM)
      +  I(CSTH*ICVM) + I(CGUI*ICVM) + 
      ICTR + I(CIDX*ICTR) + I(CUNR*ICTR) + I(CBLP*ICTR) + I(CCON*ICTR)
      +  I(CSTH*ICTR) + I(CGUI*ICTR) + 
      ICDR + I(CIDX*ICDR) + I(CUNR*ICDR) + I(CBLP*ICDR) + I(CCON*ICDR)
      +  I(CSTH*ICDR) + I(CGUI*ICDR) + 
      ICPR + I(CIDX*ICPR) + I(CUNR*ICPR) + I(CBLP*ICPR) + I(CCON*ICPR)
      + I(CSTH*ICPR) + I(CGUI*ICPR) + 
      ICJS.n + I(CIDX*ICJS.n) + I(CUNR*ICJS.n) + I(CBLP*ICJS.n) + 
      I(CCON*ICJS.n) + I(CSTH*ICJS.n) + 
      IPRI.n + I(CIDX*IPRI.n) + I(CUNR*IPRI.n) + I(CBLP*IPRI.n) + 
      I(CCON*IPRI.n) + I(CSTH*IPRI.n) + I(CGUI*IPRI.n) + 
      IDET.n + I(CIDX*IDET.n) + I(CUNR*IDET.n) + I(CBLP*IDET.n) + 
      I(CCON*IDET.n) + I(CSTH*IDET.n) + 
      IREV + I(CIDX*IREV) + I(CUNR*IREV) + I(CBLP*IREV) + I(CCON*IREV)
      +  I(CSTH*IREV) + 
      ITRI + I(CIDX*ITRI) + I(CUNR*ITRI) + I(CBLP*ITRI) + I(CCON*ITRI)
      +  I(CSTH*ITRI) + I(CGUI*ITRI),
      family = binomial)

fit2.n$linear.predictors <- (W.n %*% rowMeans(beta))[,1]
fit2.n$fitted.values <- 1/(1 + exp(-fit2.n$linear.predictors))

sample.indices <- sample(nrow(boapris.n), size=100)
boapris.ns <- boapris.n[sample.indices,]

beta.samp <- NULL
for (i in 1:86)
  beta.samp <- cbind(beta.samp,
                     eval(parse(text=paste("boapris.ns$beta.", i, ".",
                                  sep="")))) 
beta.samp <- t(beta.samp)

linpred2.n <- W.n %*% beta.samp
fitvals2.n <- list(yhat = 1/(1 + exp(-linpred2.n)))

#########
# Plots #
#########

attach("bmmp\\_Data") # bmmp functions

trellis.device()

trellis.device(postscript, file = "fig1.ps", onefile = F)

bmmp(fit1.h, samp = fitvals1.h, sparams = 6, incvar = F,
     printopts = 2)

dev.off()

bmmp(fit2.n, samp = fitvals2.n, sparams = 6, incvar = F,
     printopts = 2)

for (i in 1:39) {
  bmmp(fit1.h, samp = fitvals1.h, included = names(COUNTY[COUNTY==i]),
       sparams = 4, incvar = F, printopts = 2)}
for (i in 1:39) {
  bmmp(fit2.n, samp = fitvals2.n, included = names(COUNTY[COUNTY==i]),
       sparams = 4, incvar = F, printopts = 2)}

trellis.device(postscript, file = "fig2a.ps", onefile = F)

bmmp(fit1.h, samp = fitvals1.h, included = names(COUNTY[COUNTY==13]), 
     sparams = 4, incvar = F, printopts = 2)

dev.off()

trellis.device(postscript, file = "fig2b.ps", onefile = F)

bmmp(fit2.n, samp = fitvals2.n, included = names(COUNTY[COUNTY==13]), 
     sparams = 4, incvar = F, printopts = 2)

dev.off()

bmmp(fit1.h, samp = fitvals1.h, default = F, h = Cidx, 
     sparams = 4, incvar = F, printopts = 2)
bmmp(fit1.h, samp = fitvals1.h, default = F, h = Cunr, 
     sparams = 4, incvar = F, printopts = 2)
bmmp(fit1.h, samp = fitvals1.h, default = F, h = Cblp, 
     sparams = 4, incvar = F, printopts = 2)
bmmp(fit1.h, samp = fitvals1.h, default = F, h = Ccon, 
     sparams = 4, incvar = F, printopts = 2)

bmmp(fit2.n, samp = fitvals2.n, default = F, h = Cidx, 
     sparams = 4, incvar = F, printopts = 2)
bmmp(fit2.n, samp = fitvals2.n, default = F, h = Cunr, 
     sparams = 4, incvar = F, printopts = 2)
bmmp(fit2.n, samp = fitvals2.n, default = F, h = Cblp, 
     sparams = 4, incvar = F, printopts = 2)
bmmp(fit2.n, samp = fitvals2.n, default = F, h = Ccon, 
     sparams = 4, incvar = F, printopts = 2)

trellis.device(postscript, file = "fig3a.ps", onefile = F)

bmmp(fit1.h, samp = fitvals1.h, default = F, h = Ccon, 
     sparams = 4, incvar = F, printopts = 2)

dev.off()

trellis.device(postscript, file = "fig3b.ps", onefile = F)

bmmp(fit2.n, samp = fitvals2.n, default = F, h = Ccon, 
     sparams = 4, incvar = F, printopts = 2)

dev.off()

bmmp(fit1.h, samp = fitvals1.h, default = F, random = 1, sparams = 4,
     incvar = F) 
bmmp(fit1.h, samp = fitvals1.h, included = names(ICVS[ICVS==0]),
     sparams = 6, incvar = F) 
bmmp(fit1.h, samp = fitvals1.h, included = names(ICVS[ICVS==1]),
     sparams = 4, incvar = F) 
bmmp(fit1.h, samp = fitvals1.h, included = names(CSTH[CSTH==0]),
     sparams = 6, incvar = F) 
bmmp(fit1.h, samp = fitvals1.h, included = names(CSTH[CSTH==1]),
     sparams = 6, incvar = F) 

pprisdata.h <-
  table(Y, COUNTY)[2,]/(table(Y, COUNTY)[1,] + table(Y, COUNTY)[2,])

pprisfit1.h <- pprisfit2.n <- matrix(NA, 39, 100)

for (i in 1:39) {
  pprisfit1.h[i,] <- 
    colMeans(1/(1 + exp(-(linpred1.h[names(COUNTY[COUNTY==i]),]))))
  pprisfit2.n[i,] <- 
    colMeans(1/(1 + exp(-(linpred2.n[names(COUNTY[COUNTY==i]),]))))
}

xyplot(pprisdata.h ~ idx,
       aspect = 1,
       xlab = "Crime Rate",
       ylab = "Probability of Prison Sentence",
       panel = function(x, y, ...){
         panel.xyplot(x, y, ...)
         for (i in 1:100){
           lines(smooth.spline(x, pprisfit1.h[, i], df=3), col=4)
         }
         lines(smooth.spline(x, y, df=3), lwd=3, col=1)
       }
       )

xyplot(pprisdata.h ~ unr,
       aspect = 1,
       xlab = "Unemployment Rate",
       ylab = "Probability of Prison Sentence",
       panel = function(x, y, ...){
         panel.xyplot(x, y, ...)
         for (i in 1:100){
           lines(smooth.spline(x, pprisfit1.h[, i], df=3), col=4)
         }
         lines(smooth.spline(x, y, df=3), lwd=3, col=1)
       }
       )

trellis.device(postscript, file = "fig4a.ps", onefile = F)

xyplot(pprisdata.h ~ blp,
       aspect = 1,
       xlab = "",
       ylab = "Probability of Prison Sentence",
       panel = function(x, y, ...){
         panel.xyplot(x, y, ...)
         for (i in 1:100){
           lines(smooth.spline(x, pprisfit1.h[, i], df=3), col=4)
         }
         lines(smooth.spline(x, y, df=3), lwd=3, col=1)
       }
       )

dev.off()

xyplot(pprisdata.h ~ con,
       aspect = 1,
       xlab = "Percentage Conservative",
       ylab = "Probability of Prison Sentence",
       panel = function(x, y, ...){
         panel.xyplot(x, y, ...)
         for (i in 1:100){
           lines(smooth.spline(x, pprisfit1.h[, i], df=3), col=4)
         }
         lines(smooth.spline(x, y, df=3), lwd=3, col=1)
       }
       )

xyplot(pprisdata.h ~ idx,
       aspect = 1,
       xlab = "Crime Rate",
       ylab = "Probability of Prison Sentence",
       panel = function(x, y, ...){
         panel.xyplot(x, y, ...)
         for (i in 1:100){
           lines(smooth.spline(x, pprisfit2.n[, i], df=3), col=4)
         }
         lines(smooth.spline(x, y, df=3), lwd=3, col=1)
       }
       )

xyplot(pprisdata.h ~ unr,
       aspect = 1,
       xlab = "Unemployment Rate",
       ylab = "Probability of Prison Sentence",
       panel = function(x, y, ...){
         panel.xyplot(x, y, ...)
         for (i in 1:100){
           lines(smooth.spline(x, pprisfit2.n[, i], df=3), col=4)
         }
         lines(smooth.spline(x, y, df=3), lwd=3, col=1)
       }
       )

trellis.device(postscript, file = "fig4b.ps", onefile = F)

xyplot(pprisdata.h ~ blp,
       aspect = 1,
       xlab = "Percentage African American",
       ylab = "Probability of Prison Sentence",
       panel = function(x, y, ...){
         panel.xyplot(x, y, ...)
         for (i in 1:100){
           lines(smooth.spline(x, pprisfit2.n[, i], df=3), col=4)
         }
         lines(smooth.spline(x, y, df=3), lwd=3, col=1)
       }
       )

dev.off()

xyplot(pprisdata.h ~ con,
       aspect = 1,
       xlab = "Percentage Conservative",
       ylab = "Probability of Prison Sentence",
       panel = function(x, y, ...){
         panel.xyplot(x, y, ...)
         for (i in 1:100){
           lines(smooth.spline(x, pprisfit2.n[, i], df=3), col=4)
         }
         lines(smooth.spline(x, y, df=3), lwd=3, col=1)
       }
       )