######################################################################
###  This file is written in order to produce fuzzy wombling maps ####
###  as Figure 6 in Lu and Carlin (2005)                          ####
###  Steps involve programming in WinBUGS and Splus.              ####
######################################################################

###################################
## Step 1 WinBUGS programming    ##
###################################

model
{
   for (i in 1 : regions) {
      O[i] ~ dpois(psi[i])
      log(psi[i]) <- log(E[i]) + beta0+phi[i]
      SLDRhat[i] <- 100 * psi[i] / E[i]
   }
   phi[1:regions] ~ car.normal(adj[], weights[], num[], tau.c)
   beta0 ~ dnorm(0.0, 1.0E-5)
   tau.c ~ dgamma(0.01,0.01)
}

#Data
list(regions=87,
O =c(19,151,24,41,29,15,57,35,36,43,38,26,32,55,14,7,22,50,209,16,
42, 36,27,37,45,13,902,17,17,18,51,11,12,44,12,20,13,12,9,23,17, 28,
36,10,23,43,28,23,40,71,14,19,24,13,106,83,17,25,9,49,21, 428,6,39,28,
47,13,18,240,43,50,19,89,30,16,25,22,6,26,28,16,87,28,14, 59,63,19),

E = c(24.02630, 163.03558, 41.18794, 32.60712, 29.17479, 14.87342, 54.91725, 34.89533,
40.04383, 41.75999, 38.89972, 18.87780, 34.89533, 66.93039, 12.58520, 5.72055, 24.02630,
49.19670, 211.08817, 13.72931, 38.89972, 26.31451, 26.88657, 37.18355, 39.47177, 12.58520,
 851.21732, 24.02630, 24.59835, 19.44986, 61.20985, 12.58520, 11.44109, 41.18794, 9.15287,
 17.73369, 10.86904, 14.87342, 6.29260, 26.31451, 16.58958, 25.17040, 32.60712, 10.29698,
  20.02191, 37.18355, 28.60273, 26.31451, 41.75999, 85.23614, 10.86904, 18.87780, 26.31451,
  12.01315, 93.81696, 92.67285, 21.73808, 28.03068, 8.58082, 42.33204, 23.45424, 428.46894,
  6.86466, 29.74684, 21.73808, 44.62026, 10.29698, 18.87780, 263.71720, 44.04821, 40.61588,
   21.73808, 96.10518, 33.75122, 15.44548, 22.31013, 24.59835, 8.58082, 20.59397, 25.74246,
    18.30575, 96.10518, 23.45424, 13.72931, 53.77314, 61.78190, 25.74246),

num = c(8, 7, 7, 9, 4, 4, 7, 6, 3, 5,
8, 5, 5, 4, 6, 1, 7, 4, 6, 6,
6, 4, 4, 5, 6, 6, 7, 2, 5, 6,
5, 5, 5, 6, 2, 4, 4, 2, 3, 6,
4, 5, 5, 4, 5, 5, 5, 7, 6, 5,
8, 5, 5, 4, 6, 7, 5, 5, 5, 6,
6, 4, 2, 6, 9, 7, 3, 4, 5, 6,
7, 6, 9, 6, 6, 6, 6, 4, 3, 5,
6, 4, 5, 4, 4, 7, 6
),
weights=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1),

adj = c(
69, 58, 48, 33, 31, 18, 11, 9,
86, 82, 71, 62, 30, 27, 13,
80, 56, 54, 44, 29, 15, 14,
68, 57, 45, 39, 36, 31, 29, 15, 11,
73, 71, 49, 48,
78, 76, 75, 37,
83, 81, 52, 46, 40, 22, 8,
83, 65, 64, 52, 17, 7,
69, 58, 1,
86, 72, 70, 43, 27,
80, 77, 49, 31, 29, 18, 4, 1,
87, 76, 65, 37, 34,
82, 58, 33, 30, 2,
84, 56, 54, 3,
60, 57, 44, 29, 4, 3,
38,
83, 64, 53, 51, 46, 32, 8,
49, 48, 11, 1,
82, 70, 66, 62, 27, 25,
74, 66, 55, 50, 25, 24,
77, 75, 73, 61, 56, 26,
81, 46, 24, 7,
85, 55, 50, 28,
81, 74, 50, 22, 20,
79, 74, 66, 55, 20, 19,
84, 78, 75, 61, 56, 21,
86, 71, 70, 62, 19, 10, 2,
85, 23,
80, 15, 11, 4, 3,
71, 58, 48, 33, 13, 2,
69, 36, 11, 4, 1,
83, 53, 51, 46, 17,
58, 48, 30, 13, 1,
76, 73, 65, 61, 47, 12,
68, 45,
69, 39, 31, 4,
87, 76, 12, 6,
69, 16,
68, 36, 4,
81, 72, 70, 66, 52, 7,
87, 59, 51, 42,
87, 64, 59, 51, 41,
86, 72, 65, 47, 10,
60, 54, 15, 3,
68, 60, 57, 35, 4,
83, 32, 22, 17, 7,
86, 73, 65, 43, 34,
71, 49, 33, 30, 18, 5, 1,
77, 73, 48, 18, 11, 5,
74, 55, 24, 23, 20,
67, 64, 59, 53, 42, 41, 32, 17,
72, 65, 40, 8, 7,
67, 59, 51, 32, 17,
60, 44, 14, 3,
85, 79, 50, 25, 23, 20,
84, 80, 77, 26, 21, 14, 3,
63, 60, 45, 15, 4,
33, 30, 13, 9, 1,
67, 53, 51, 42, 41,
63, 57, 54, 45, 44, 15,
76, 75, 73, 34, 26, 21,
82, 27, 19, 2,
60, 57,
87, 65, 51, 42, 17, 8,
87, 72, 64, 52, 47, 43, 34, 12, 8,
81, 74, 70, 40, 25, 20, 19,
59, 53, 51,
45, 39, 35, 4,
38, 36, 31, 9, 1,
72, 66, 40, 27, 19, 10,
86, 73, 48, 30, 27, 5, 2,
70, 65, 52, 43, 40, 10,
86, 77, 71, 61, 49, 47, 34, 21, 5,
81, 66, 50, 25, 24, 20,
78, 76, 61, 26, 21, 6,
75, 61, 37, 34, 12, 6,
80, 73, 56, 49, 21, 11,
84, 75, 26, 6,
85, 55, 25,
77, 56, 29, 11, 3,
74, 66, 40, 24, 22, 7,
62, 19, 13, 2,
46, 32, 17, 8, 7,
78, 56, 26, 14,
79, 55, 28, 23,
73, 71, 47, 43, 27, 10, 2,
65, 64, 42, 41, 37, 12
))


#Inits
list(tau.c = 1, beta0 = 0, phi=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
))

###########################################################
## Note: save the coda of SLDRhat for Splus programming ###
## and saved as "sldrhat.txt"                           ###
###########################################################


######################################
##  Step 2  Splus programming      ###
######################################

############################################################
## Read in necessary data such as adjacency information   ##
############################################################


fixed.map.ng.num<-c(8,6,7,9,4,4,7,6,3,5,
8,5,4,4,6,1,5,4,6,4,
5,4,4,4,5,5,6,2,5,5,
5,3,4,6,2,4,4,2,3,6,
3,4,5,4,5,4,5,7,6,4,
5,5,3,4,6,7,5,4,3,6,
5,4,2,6,9,6,2,4,5,6,
6,6,9,4,5,6,6,4,3,5,
6,4,4,4,4,6,6)

fixed.map.adj<-c(
69, 58, 48, 33, 31, 18, 11, 9,
82, 71, 62, 30, 27, 13,
80, 56, 54, 44, 29, 15, 14,
68, 57, 45, 39, 36, 31, 29, 15, 11,
73, 71, 49, 48,
78, 76, 75, 37,
83, 81, 52, 46, 40, 22, 8,
83, 65, 64, 52, 17, 7,
69, 58, 1,
86, 72, 70, 43, 27,
80, 77, 49, 31, 29, 18, 4, 1,
87, 76, 65, 37, 34,
82, 58,  30, 2,
84, 56, 54, 3,
60, 57, 44, 29, 4, 3,
38,
83, 64,  51,  32, 8,
49, 48, 11, 1,
82, 70, 66, 62, 27, 25,
74,  55, 50, 25,
77,  73, 61, 56, 26,
81, 46, 24, 7,
85, 55, 50, 28,
81, 74, 50, 22,
79,  66, 55, 20, 19,
84, 78, 75,  56, 21,
86,  70, 62, 19, 10, 2,
85, 23,
80, 15, 11, 4, 3,
71,  48, 33, 13, 2,
69, 36, 11, 4, 1,
53,  46, 17,
58, 48, 30,  1,
76, 73, 65, 61, 47, 12,
68, 45,
69, 39, 31, 4,
87, 76, 12, 6,
69, 16,
68, 36, 4,
81, 72, 70, 66, 52, 7,
87, 59,  42,
87, 64, 51, 41,
86, 72, 65, 47, 10,
60, 54, 15, 3,
68, 60, 57, 35, 4,
83, 32, 22,  7,
86, 73, 65, 43, 34,
71, 49, 33, 30, 18, 5, 1,
77, 73, 48, 18, 11, 5,
 55, 24, 23, 20,
64, 59, 53, 42,   17,
72, 65, 40, 8, 7,
67,  51, 32,
60, 44, 14, 3,
85, 79, 50, 25, 23, 20,
84, 80, 77, 26, 21, 14, 3,
63, 60, 45, 15, 4,
33,  13, 9, 1,
67,  51, 41,
63, 57, 54, 45, 44, 15,
76, 75, 73, 34, 21,
82, 27, 19, 2,
60, 57,
87, 65, 51, 42, 17, 8,
87, 72, 64, 52, 47, 43, 34, 12, 8,
81, 74, 70, 40, 25, 19,
59, 53,
45, 39, 35, 4,
38, 36, 31, 9, 1,
72, 66, 40, 27, 19, 10,
86, 73, 48, 30,  5, 2,
70, 65, 52, 43, 40, 10,
86, 77, 71, 61, 49, 47, 34, 21, 5,
81, 66, 24, 20,
78, 76, 61, 26,  6,
75, 61, 37, 34, 12, 6,
80, 73, 56, 49, 21, 11,
84, 75, 26, 6,
85, 55, 25,
77, 56, 29, 11, 3,
74, 66, 40, 24, 22, 7,
62, 19, 13, 2,
46, 17, 8, 7,
78, 56, 26, 14,
79, 55, 28, 23,
73, 71, 47, 43, 27, 10,
65, 64, 42, 41, 37, 12
)

ind1<-c(1,9,15,22,31,35,39,46,52,55,60,68,73,77,81,87,88,93,97,103,107,112,116,
120,124,129,134,140,142,147,152,157,160,164,170,172,176,180,182,185,191
,194,198,203,207,212,216,221,228,234,238,243,248,251,255,261,268,273,27
7,280,286,291,295,297,303,312,318,320,324,329,335,341,347,356,360,365,3
71,377,381,384,389,395,399,403,407,411,417)

ind2<-c(8,14,21,30,34,38,45,51,54,59,67,72,76,80,86,87,92,96,102,106,111,115,11
9,123,128,133,139,141,146,151,156,159,163,169,171,175,179,181,184,190,1
93,197,202,206,211,215,220,227,233,237,242,247,250,254,260,267,272,276,
279,285,290,294,296,302,311,317,319,323,328,334,340,346,355,359,364,370
,376,380,383,388,394,398,402,406,410,416,422)



#############################################
# First compute Delta=|SLDRhati-SLDRhatj|   #
#############################################

temp<-read.table("sldrhat.txt")[,2]
sldrhat<-matrix(0, ncol=87, nrow=2000) ## if 2000 is the number of production iteration in WinBUGS
for(i in 1:87){
    from<-(i-1)*2000+1
    to<-i*2000
    sldrhat[,i]<-temp[from:to]
}

delta<-matrix(0, ncol=sum(fixed.map.ng.num)/2,nrow=2000)
k<-0
for( i in 1:87){
    for(j in ind1[i]:ind2[i]){
        if(fixed.map.adj[j]>i){
            k<-k+1
            delta[,k]<-abs(sldrhat[,i]-sldrhat[,fixed.map.adj[j]])
            }}}

#############################################
# Then prepare functions and data needed for#
# draw the wombling maps in Splus           #
#############################################

library(maps)
ob.mn <- map('county','minnesota',fill=F,add=T,plot=F)

index <- seq(1,length(ob.mn$x))
index.na <-index[ob.mn$x=="NA"]
from.ind <- rep(0,length(index.na)+1)
to.ind <- rep(0,length(index.na)+1)
from.ind[1] <- 1
to.ind[1] <- 2
from.ind[length(from.ind)] <- 1562
to.ind[length(from.ind)] <- 1563

for ( i in 2:length(index.na)){
  from.ind[i] <-index[(index.na[i-1]+1)]
  to.ind[i] <- index[(index.na[i]-1)]
}

mn.margin.ind<-c(162,254,174,165,257,171,86,173,
256,255,54,234,233,72,76,275,274,273,94,128,127,271,
270,277,279,278,138,137,141,140,118,119,212,211,210,
124,123,113,112,200,202,152,151,150,219,221,253,252,
236,185,184,280,168,167,38,37,269,268,276,78,224,223,240,239,196,195,163)

from.v1<-from.ind[-c(235,mn.margin.ind)]
to.v1<-to.ind[-c(235,mn.margin.ind)]

map.lines.order<-c(5,7,6,3,4,2,8,1,13,11,9,14,12,10,15,19,20,18,17,21,16,24,23,29,28,30,27,26,25,22,
33,31,32,34,37,35,38,36,42,44,40,43,41,45,39,47,48,50,46,49,52,51,57,53,56,54,55,62,58,59,60,63,61,
65,67,64,66,68,70,71,69,74,72,73,75,77,78,76,79,81,82,83,80,84,85,91,86,89,90,87,88,94,93,95,92,97,
96,99,100,98,103,101,102,107,104,105,106,110,109,108,113,111,112,117,114,115,116,120,119,118,121,
122,125,124,123,126,127,129,128,130,131,135,132,134,136,133,138,137,139,140,141,142,143,144,
146,149,145,148,147,151,152,150,153,155,154,157,156,159,158,160,161,163,164,162,165,167,168,166,
169,170,171,172,173,175,176,174,177,178,179,180,181,182,185,184,183,187,186,188,189,190,193,192,191,
194,196,195,197,198,201,199,200,202,204,203,205,206,207,208,209,210,211,212)

ordered.from<-ordered.to<-rep(0,212)
ordered.from[map.lines.order]<-from.v1
ordered.to[map.lines.order]<-to.v1

use.delta.order<-c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,
26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,52,53,55,
56,57,58,59,60,61,62,63,64,65,68,69,70,71,72,73,74,75,77,78,79,81,82,83,84,87,88,89,
90,91,93,94,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,116,
117,118,119,120,121,122,124,125,126,129,130,131,132,134,135,136,140,142,147,148,149,
152,153,157,158,160,161,164,165,166,167,168,170,171,172,173,176,177,180,182,185,186,
187,188,189,191,192,193,194,195,196,198,199,200,201,203,204,207,208,209,212,216,217,
218,221,222,228,229,234,238,239,240,243,244,248,251,255,256,261,262,263,268,269,277,
280,286,287,288,291,297,298,303,304,312,313,314,329,335,336,347,348,356,360,361,371,377,381)

draw.var<-function(x,divider,lgd=T,t){
    if(length(x)!=211){x<-x[use.delta.order]}
    sq_seq(1,211,1)

    ln1<-sq[as.numeric(x<=divider[2])==1]
    ln2<-sq[as.numeric (x>divider[2]& x<=divider[3])==1]
    ln3<-sq[as.numeric (x>divider[3]& x<=divider[4])==1]
    ln4<-sq[as.numeric (x>divider[4]& x<=divider[5])==1]

    color<-0*1:211
    color[ln1] <- 13
    color[ln2] <- 0
    color[ln3] <- 15
    color[ln4] <- 15

    od<-function(i,c,wide){lines(ob.mn$x[ordered.from[i]:ordered.to[i]],ob.mn$y[ordered.from[i]:ordered.to[i]],col=c,lwd=wide)}
    draw.line<-function(ln,c,wide){
        for ( i in 1:length(ln)){
            if(ln[i]<187){od(ln[i],c,wide)}
            if(ln[i]==187){od(ln[i],c,wide);od(ln[i]+1,c,wide)}
            if(ln[i]>187){od(ln[i]+1,c,wide)}
    }
    }
    map('county','minnesota',fill=T,add=F,col=12)
    map('county','minnesota',fill=F,add=T,col=2,lwd=0.8)
    draw.line(ln1,0,5)
    draw.line(ln2,15,5)
    draw.line(ln3,3,5)
    draw.line(ln4,1,5)

    if(lgd==T){legend(-92,46, c(paste(round(divider[1],2),"-",round(divider[2],2)),paste(round(divider[2],2),"-",round(divider[3],2)),paste(round(divider[3],2),"-",round(divider[4],2)), paste(round(divider[4],2),"-",round(divider[5],2))),fill=c(0,15,3,1),bty="n")}
    title(t)
}


#################################################
### now compute phat sd(phat) with different   ##
### cutoff point c=5, 15, and 30               ##
#################################################

c.func5<-function(x){
    c.ind<-as.numeric(x>5)
    return(c.ind)
}

p.y.5<-apply(apply(delta,2,c.func5)/2000,2,sum)[use.delta.order]
sd.p.y.5<-sqrt(p.y.5*(1-p.y.5)/2000)

c.func15<-function(x){
    c.ind<-as.numeric(x>15)
    return(c.ind)
}

p.y.15<-apply(apply(delta,2,c.func15)/2000,2,sum)[use.delta.order]
sd.p.y.15<-sqrt(p.y.15*(1-p.y.15)/2000)


c.func30<-function(x){
    c.ind<-as.numeric(x>30)
    return(c.ind)
}

p.y.30<-apply(apply(delta,2,c.func30)/2000,2,sum)[use.delta.order]
sd.p.y.30<-sqrt(p.y.30*(1-p.y.30)/2000)


draw.var(p.y.5,c(0,0.4,0.6,0.8,1),t="phat, c=5")
draw.var(p.y.15,c(0,0.4,0.6,0.8,1),t="phat, c=15")
draw.var(p.y.30,c(0,0.4,0.6,0.8,1),t="phat, c=30")


draw.var(sd.y.5,c(0,0.006,0.008,0.01,0.012),t="sd(phat), c=5")
draw.var(sd.y.15,c(0,0.006,0.008,0.01,0.012),t="sd(phat), c=15")
draw.var(sd.y.30,c(0,0.006,0.008,0.01,0.012),t="sd(phat), c=30")