library("arm")
   setwd("C:/Documents and Settings/zhihengxu/My Documents/Ph.D. research/Computations/AD")
   aapc <- read.table ("iowa.txt", header=TRUE)
   adj<-read.table("IA_counties_adj.txt")[,1]
   num<-read.table("IA_counties_num.txt")[,1]
   Nneigh<-588
   A <- 10
   P<-6
   Co<-15
   C<-99
   county <- aapc[,1]
   age <- aapc[,2]
   period <- aapc[,3]
   cohort <- A+period-age
   pyr <- aapc[,4]
   cases <- aapc[,5]
   rate<-100000*cases/pyr

#create adjacent matrix
adj_max<-mat.or.vec(C,C)
begin<-1
end<-0
for (i in 1:C){
	end<-end+num[i]
	for (j in begin:end){
		adj_max[i,adj[j]]<-1
	}
	begin<-begin+num[i]
}


vec10<-c(1:10)*0.01
vec6<-c(1:6)*0.01
vec99<-c(1:99)*0.01
vec15<-c(1:15)*0.01
 
# area and age
   data <- list ("Nneigh","adj","num", "C","county","age", "A","cases", "pyr")
   inits <- function() {list (tauc=1, taua=1, r=0.07, p2=-0.0885, q=0.000029, alpha=vec10,u=vec99,v=vec99)}
   parameters <- c("tauc", "u","alphac","taua","v", "r", "p2", "q", "Dp")
   aapc.sim <- bugs (data, inits, parameters, "aa_pred.bug", n.chains=1, n.iter=1000, debug=TRUE)

# area and age,period
   data <- list ("Nneigh","adj","num", "C","county","age", "A","period","P","cases", "pyr")
   inits <- function() {list (tauc=1, taua=1,taua2=10, s=5,c=-10, u=vec99,alphac=vec10,v=vec99,beta=vec6)}
   parameters <- c("tauc", "u","alphac","taua2", "taua", "c", "s","taup","betac")
   aapc.sim <- bugs (data, inits, parameters, "aap.bug", n.chains=1, n.iter=100, debug=TRUE)

# area, age, and cohort
   data <- list ( "Nneigh","adj","num","C","county","age", "A", "cohort","Co","cases", "pyr")
   inits <- function() {list ( taua=1, tauco=1,r=0.07, p2=-0.0885, q=0.000029,alpha=vec10, u=vec99,v=vec99,gamma=vec15)}
   parameters <- c( "u","v","alphac","r", "p2", "q","taua", "gammac","tauco")
   aapc.sim <- bugs (data, inits, parameters, "aac.bug", n.chains=1, n.iter=1000, debug=TRUE)

# area and age,period and area-period interaction
   data <- list ("Nneigh","adj","num","C","county","age", "A", "period", "P", "cases", "pyr")
   inits <- function() {list ( taua=1, taup=1,r=0.07, p2=-0.0885, q=0.000029, alpha=vec10, beta=vec6,
			u=vec99,v=vec99)}
   parameters <- c("tauc", "u","v","alphac","r", "p2", "q","taua", "betac", "taup", "phi_area","phi_period")
   aapc.sim <- bugs (data, inits, parameters, "aap+ap.bug", n.chains=1, n.iter=10000, debug=TRUE)

# area, age,period and cohort
   data <- list ( "Nneigh","adj","num","C","county","age", "A", "period", "P","cohort","Co","cases", "pyr")
   inits <- function() {list ( taua=1, taup=1,tauco=1,r=0.07, p2=-0.0885, q=0.000029, alpha=vec10, beta=vec6,u=0,v=vec99,gamma=vec15)}
   parameters <- c( "u","v","alphac","r", "p2", "q","taua", "betac", "taup","gammac","tauco", "phi","phi2")
   aapc.sim <- bugs (data, inits, parameters, "aapc.bug", n.chains=1, n.iter=10000, debug=TRUE)


# AAPC + AP
   data <- list ( "Nneigh","adj","num","C","county","age", "A", "period", "P","cohort","Co","cases", "pyr")
   inits <- function() {list ( taua=1, taup=1,tauco=1,r=0.07, p2=-0.0885, q=0.000029, alpha=vec10, beta=vec6,u=0,v=vec99,gamma=vec15)}
   parameters <- c( "u","v","alphac","r", "p2", "q","taua", "betac", "taup","gammac","tauco", "phi_area","phi_period")
   aapc.sim <- bugs (data, inits, parameters, "aapc+ap.bug", n.chains=1, n.iter=1000, debug=TRUE)

# AAPC + AC
   data <- list ( "Nneigh","adj","num","C","county","age", "A", "period", "P","cohort","Co","cases", "pyr")
   inits <- function() {list ( taua=1, tau2=1,taup=1,tauco=1,cons=-10,s=5, alpha=vec10, beta=vec6,u=0,v=vec99,gamma=vec15)}
   parameters <- c( "u","v","alphac","s","cons","taua", "betac", "taup","gammac","tauco", "phi_area","phi_cohort")
   aapc.sim <- bugs (data, inits, parameters, "aapc+ac.bug", n.chains=1, n.iter=100, debug=TRUE, save.history=TRUE)

# AAPC + AC2
   data <- list ( "Nneigh","adj","num","C","county","age", "A", "period", "P","cohort","Co","cases", "pyr")
   inits <- function() {list ( taua=1, tau2=1,taup=1,tauco=1,alpha=vec10, beta=vec6,u=0,v=vec99,gamma=vec15)}
   parameters <- c( "u","v","alphac","taua", "betac", "taup","gammac","tauco", "phi_area","phi_cohort")
   aapc.sim <- bugs (data, inits, parameters, "aapc+ac2.bug", n.chains=1, n.iter=1000, debug=TRUE, save.history=TRUE)


# AAPC + AP+ AC
   data <- list ( "Nneigh","adj","num","C","county","age", "A", "period", "P","cohort","Co","cases", "pyr")
   inits <- function() {list ( taua=1, taup=1,tauco=1, s=5, cons=0,alpha=vec10, beta=vec6,u=0,v=vec99,gamma=vec15)}
   parameters <- c( "u","v","alpha","alphac","taua", "s", "cons","betac", "taup","gammac","tauco", "phi_area","phi_cohort")
   aapc.sim <- bugs (data, inits, parameters, "aapc+apac.bug", n.chains=1, n.iter=1000, debug=TRUE, save.history=TRUE,codaPkg=TRUE)

# AAPC (car.normal for period and cohort)
   data <- list ( "Nneigh","adj","num","C","county","age", "A", "period", "P","cohort","Co","cases", "pyr")
   inits <- function() {list ( taua=1, taup=1, tauc=1, tauco=1,s=5, cons=0,alpha=vec10, beta=vec6,u=0,v=vec99,gamma=vec15)}
   parameters <- c( "u","v","alpha","alphac", "s", "cons","betac","gammac")
   aapc.sim <- bugs (data, inits, parameters, "aapc2.bug", n.chains=1, n.iter=100, debug=TRUE, save.history=TRUE)


# area and age,cohort and area-cohort interaction
   data <- list ("Nneigh","adj","num","C","county","age", "A", "cohort", "Co", "cases", "pyr")
   inits <- function() {list ( taua=1, tauco=1,r=0.07, p2=-0.0885, q=0.000029, alpha=vec10, gamma=vec15,
			u=0,v=vec99)}
   parameters <- c("tauc", "u","v","alphac","r", "p2", "q","taua", "gammac", "tauco", "phi","phi_period")
   aapc.sim <- bugs (data, inits, parameters, "aacac.bug", n.chains=1, n.iter=1000, debug=TRUE)

# area and age,cohort and area-cohort interaction
   data <- list ("Nneigh","adj","num","C","county","age", "A", "cohort", "Co", "cases", "pyr")
   inits <- function() {list ( taua=1, tauco=1, alpha=vec10, gamma=vec15,
			u=0,v=vec99)}
   parameters <- c("tauc", "u","v","alphac","taua", "gammac", "tauco", "phi","phi_period")
   aapc.sim <- bugs (data, inits, parameters, "aapc_aacACold.bug", n.chains=1, n.iter=1000, debug=TRUE)



# area, age,period and area-period interaction
   data <- list ( "adj_max","Nneigh","adj","num","C","county","age", "A", "period", "P", "cases", "pyr")
   inits <- function() {list (tauc=1, taua=1, taup=1,tauphi=1,s=5, alpha=mat.or.vec(10,1), beta=mat.or.vec(6,1),
			u=0,v=mat.or.vec(99,1), phi=mat.or.vec(99,6)+0.001)}
   parameters <- c("tauc", "u","v","alphac","s","taua", "betac", "taup", "phi","tauphi")
   aapc.sim <- bugs (data, inits, parameters, "aapc_inter.bug", n.chains=1, n.iter=1000, debug=TRUE)


# area, age,period and area-period interaction-Condgon formula
   data <- list ("Nneigh","adj","num","C","county","age", "A", "period", "P", "cases", "pyr")
   inits <- function() {list (tauc=1, taua=1, taup=1,s=5, alpha=vec10, beta=vec6,
			u=0,v=vec99, phi=mat99by6)}
   parameters <- c("tauc", "u","v","alphac","s","taua", "betac", "taup", "phi","phi_period")
   aapc.sim <- bugs (data, inits, parameters, "aapc_condgon.bug", n.chains=1, n.iter=1000, debug=TRUE)