###########################################################
## R code used to produce plots in Ma, Virnig and Carlin ##
##		08/15/2005 Haijun Ma			 ##
###########################################################

#################################################################################
## Need to install R package "maps" first. Download it from www.r-project.org	#
## In R console, first call in R functions wrote for producing MN county maps.  #
## The functions are saved at the current directory in the file "MVCmnPlot.txt".#
#################################################################################


source("MVCmnPlot.txt")




###################################################
###   Figure 1. MN Breast Cancer raw Data       ###
###################################################

# save the raw data in the current directory in the file DATA.txt.
source("DATA.txt")

postscript(file = "SE_bcA_rSMR.eps", onefile = T,paper="special", width=5, height=6, horizontal=F)
par(mar=c(1,2,3,2)+0.1)
mn.Poly(expression(paste("(a) Raw SLDR", sep="")),SMR , 5)
dev.off()

postscript(file = "SE_bcA_rD.eps", onefile = T,paper="special", width=5, height=6, horizontal=F)
par(mar=c(1,2,3,2)+0.1)
mn.Edge.crisp(expression(paste("(b) Raw ", D[ij],"=", abs(SLDR[i]-SLDR[j]),sep="")),D , T, c(0,.6,.8,.9,1))
dev.off()

postscript(file = "SE_bcA_rmamg.eps", onefile = T,paper="special", width=5, height=6, horizontal=F)
par(mar=c(1,2,3,2)+0.1)
mn.Poly(expression(paste("(c) Nonscreening rate, %", sep="")),100-origX , 5)
dev.off()

###################################################
###   Figure 2.MN Breast Cancer Global smooting ###
###################################################

# WinBUGS output saved in the current directory as .txt files.

data<-read.table("AbcEta_areaIID.txt", skip=1, header=F)
RR<-data[,2]
data<-read.table("AbcDelta_areaIID.txt", skip=1, header=F)
delta<-data[,2]
data<-read.table("AbcDeltamy_areaIID.txt", skip=1, header=F)
deltamy<-data[,2]

postscript(file = "bcA_eRR_iid.eps", onefile = T,paper="special", width=5, height=6, horizontal=F)
par(mar=c(1,2,3,2)+0.1)
mn.Poly(expression(paste("(a) E( ",eta," | ",bolditalic(y)," )")),RR , 5)
dev.off()

postscript(file = "bcA_eD_iid.eps", onefile = T,paper="special", width=5, height=6, horizontal=F)
par(mar=c(1,2,3,2)+0.1)
mn.Edge.crisp(expression(paste("(b) E( |",eta[i],"-",eta[j],"| | ",bolditalic(y)," )",sep="")),delta,T,c(0,.6,.8,.9,1))
dev.off()

postscript(file = "bcA_eDmy_iid.eps", onefile = T,paper="special", width=5, height=6, horizontal=F)
par(mar=c(1,2,3,2)+0.1)
mn.Edge.crisp(expression(paste("(c) | E(",eta[i],"-",eta[j],")| ",bolditalic(y)," |",sep="")),abs(deltamy),T,c(0,.6,.8,.9,1))
dev.off()

############################################
##  Figure 3 (a). explain area adjacency  ##
############################################
par(mfrow=c(1,1))
plot(xy.coord[,1],xy.coord[,2],xlab=NA, ylab=NA, axes=F,type="n", sub="(a) Area Adjacency")
polygon(xy.coord[,1],xy.coord[,2],col=c("white", "white","grey","white","white",
"grey","gray30","grey","white","white",
"grey","white","white","white","white","white"))

# lable the grids
legend.x<-rep(c(1.5,2.5,4,6),4)
legend.y<-c(rep(1.5,4),rep(3,4),rep(4.5,4),rep(5.5,4))
for(i in 1:length(legend.x)){
legend(legend.x[i],legend.y[i],legend=as.character(i),bty="n",xjust=0.7,yjust=0.5, cex=2)
}
dev.print(device=postscript, "AreaAdj.eps",onefile=FALSE, horizontal=FALSE, paper="special")


##########################################
##  Figure 3 (b). explain edge adjacency #
##########################################

par(mfrow=c(1,1))
plot(xy.coord[,1],xy.coord[,2],xlab=NA, ylab=NA, axes=F,type="n", sub="(b) Edge Adjacency")
polygon(xy.coord[,1],xy.coord[,2])
# lable the grids
legend.x<-rep(c(1.5,2.5,4,6),4)
legend.y<-c(rep(1.5,4),rep(3,4),rep(4.5,4),rep(5.5,4))
for(i in 1:length(legend.x)){
legend(legend.x[i],legend.y[i],legend=as.character(i),bty="n",xjust=0.7,yjust=0.5,cex=2)
}
# highlight edge(6,11)
segments(3,4,5,4, col="black",lwd=7)
#end.coords give the starting and ending coords for each segment. 4 values per segment
end.coords<-c(2,4,3,4,3,4,3,2,5,4,7,4,5,4,5,2,3,5,3,4,5,5,5,4)
# calculate number of segments
n.edge<-length(end.coords)/4
end.y<-end.coords[4*1:6]
end.x<-end.coords[4*1:6-1]
st.y<-end.coords[4*1:6-2]
st.x<-end.coords[4*1:6-3]
segments(st.x,st.y,end.x,end.y, col="black", lwd=7, lty="dashed")
dev.print(device=postscript, "EdgeAdj.eps",onefile=FALSE, horizontal=FALSE, paper="special")

####################################################
###   Figure 4. MN Breast Cancer Local Smoothing ###
####################################################

# WinBUGS output saved in the current directory as .txt files.

data<-read.table("AbcEta_areaCAR.txt", skip=1, header=F)
RR<-data[,2]
data<-read.table("AbcDelta_areaCAR.txt", skip=1, header=F)
delta<-data[,2]
data<-read.table("AbcDeltamy_areaCAR.txt", skip=1, header=F)
deltamy<-data[,2]


postscript(file = "bcA_eRR_CAR.eps", onefile = T,paper="special", width=5, height=6, horizontal=F)
par(mar=c(1,2,3,2)+0.1)
mn.Poly(expression(paste("(a) E( ",eta," | ",bolditalic(y)," )")),RR , 5)
dev.off()

postscript(file = "bcA_eD_CAR.eps", onefile = T,paper="special", width=5, height=6, horizontal=F)
par(mar=c(1,2,3,2)+0.1)
mn.Edge.crisp(expression(paste("(b) E( |",eta[i],"-",eta[j],"| | ",bolditalic(y)," )",sep="")),delta,T,c(0,.6,.8,.9,1))
dev.off()

postscript(file = "bcA_eDmy_CAR.eps", onefile = T,paper="special", width=5, height=6, horizontal=F)
par(mar=c(1,2,3,2)+0.1)
mn.Edge.crisp(expression(paste("(c) | E(",eta[i],"-",eta[j],")| ",bolditalic(y)," | ",sep="")),abs(deltamy),T,c(0,.6,.8,.9,1))
dev.off()



##########################################################
###   Figure 5. MN Breast Cancer Local Edge Smoothing  ###
##########################################################

# function used to add in the random effects for the two "island" edges
insert<-function(phiE,phi){
temp<-seq(0,211)
temp[1:78]<-phiE[1:78]
temp[79]<-phi[1]
temp[80:142]<-phiE[79:141]
temp[143]<-phi[2]
temp[144:211]<-phiE[142:209]
return(temp)
}

# WinBUGS output saved in the current directory as .txt files.

data<-read.table("AbcEta_ED.txt", skip=1, header=F)
eta<-data[,2]
data<-read.table("AbcPsi_ED.txt", skip=1, header=F)
psi<-data[,2]
data<-read.table("AbcPhi_ED.txt", skip=1, header=F)
phi<-data[,2]
psimean<-insert(psi,phi)

index<-scan("AbcIndex_ED.txt",sep=",")
ind1<-index[1:211]
ind2<-index[212:422]
ind3<-index[423:633]


postscript(file = "bcA_Eeta_DE.eps", onefile = T,paper="special", width=5, height=6, horizontal=F)
par(mar=c(1,2,3,2)+0.1)
mn.Edge.crisp(expression(paste("(a) E( ",delta[ij], " | ",bolditalic(y), ")",sep="")),eta,T,c(0,.6,.8,.9,1))
dev.off()


postscript(file = "bcA_Epsi_DE.eps", onefile = T,paper="special", width=5, height=6, horizontal=F)
par(mar=c(1,2,3,2)+0.1)
mn.Edge.crisp(expression(paste("(b) E( ",psi[ij], " | ", bolditalic(y)," )",sep="")),psimean,T,c(0,.6,.8,.9,1))
dev.off()



postscript(file = "bcA_PrD1_DE.eps", onefile = T,paper="special", width=5, height=6, horizontal=F)
par(mar=c(1,2,3,2)+0.1)
mn.Edge.crisp(expression(paste("(c) Pr( ", delta[ij], " > -0.161 | ", bolditalic(y)," )",sep="")),ind1,F,c(0,.6,.8,.9,1))
dev.off()

postscript(file = "bcA_PrD2_DE.eps", onefile = T,paper="special", width=5, height=6, horizontal=F)
par(mar=c(1,2,3,2)+0.1)
mn.Edge.crisp(expression(paste("(d) Pr( ", delta[ij], " > -1.39 | ", bolditalic(y)," )",sep="")),ind2,F,c(0,.6,.8,.9,1))
dev.off()


postscript(file = "bcA_PrD3_DE.eps", onefile = T,paper="special", width=5, height=6, horizontal=F)
par(mar=c(1,2,3,2)+0.1)
mn.Edge.crisp(expression(paste("(e) Pr( ", delta[ij], " > -1.20 | ", bolditalic(y)," )",sep="")),ind3,F,c(0,.6,.8,.9,1))
dev.off()