###################Read data(y, x1, x2, Time, pop_tn, rate)##################
library(dplyr) 
library(readr)
setwd("data")
############ x1 - unemployment rate
unemploy <- matrix(0, nrow = 83, ncol = 14)
for (i in 1:14) {
  year <- formatC(c(90:96, 98:99, 00:04), width=2, flag = 0)
  unemploy[, i] <- read_csv(sprintf("laucnty%s.csv", year[i])) %>% 
    filter(state == 26) %>% select(Unemployment_rate) %>% 
    mutate(Unemployment_rate/100) %>% pull
}
xx1 <- as.vector(unemploy)


############ x2 - rurality
rurality <- read.table("injury.txt", header=T) %>% 
  filter(YEAR == 1990) %>% select(RURALITY) %>% pull
xx2 <- rep(rurality, 14)


############ Time = year_j - 1997
Time <- c(-7:-1, 1:7) %>% rep(each=83)

############ y - teenage car crash counts
yy <- read_csv("mi_GDL1617.csv") %>% 
  filter(YEAR %in% c(1990:1996, 1998:2004)) %>%
  select(count_tn) %>% pull %>% 
  matrix(83, 14, byrow = T) %>% as.vector

############ n - pop_tn
pop_tn <- read_csv("mi_GDL1617.csv") %>% 
  filter(YEAR %in% c(1990:1996, 1998:2004)) %>%
  select(pop_tn) %>% pull %>% 
  matrix(83, 14, byrow = T) %>% as.vector

############ teen driver&rate of teen driver(state level)
pop_tn_state <- read_csv("mi_GDL1617.csv") %>% 
  filter(YEAR %in% c(1990:1996, 1998:2004)) %>%
  group_by(YEAR) %>% summarise(state=sum(pop_tn)) %>% pull(state)
pop_tn_state_plot <- read_csv("mi_GDL1617.csv") %>% 
  filter(YEAR %in% c(1990:2004)) %>%
  group_by(YEAR) %>% summarise(state=sum(pop_tn)) %>% pull(state)
teen_d <- read.table("mi_teen_pop.txt", header = T) %>%
  filter(year %in% c(1990:1996, 1998:2004)) %>% select(under19) %>% pull
teen_d_plot <- read.table("mi_teen_pop.txt", header = T) %>%
  filter(year %in% c(1990:2004))
rate <- teen_d/pop_tn_state
#########################Finish reading data here############################


##################Get Arguments for CAR prior################################
library(maptools); library(glmmBUGS); library(maps)
############ transform the dataset from maps into a SpatialPolygon object
mi.county <- maps::map("county", "michigan", fill = T, plot = F)
mi.poly.county <- maptools::map2SpatialPolygons(mi.county,
                                                IDs=mi.county$names)
car.values <- glmmBUGS::addSpatial(mi.poly.county)
adj <- car.values$adj
num <- car.values$num
sumNumNeigh <- sum(num)

############ transfer function, from bugs to stan (get adj matrix)
mungeCARdata4stan <- function(adjBUGS,numBUGS) {
  N = length(numBUGS);
  nn = numBUGS;
  N_edges = length(adjBUGS) / 2;
  node1 = vector(mode="numeric", length=N_edges);
  node2 = vector(mode="numeric", length=N_edges);
  iAdj = 0;
  iEdge = 0;
  for (i in 1:N) {
    for (j in 1:nn[i]) {
      iAdj = iAdj + 1;
      if (i < adjBUGS[iAdj]) {
        iEdge = iEdge + 1;
        node1[iEdge] = i;
        node2[iEdge] = adjBUGS[iAdj];
      }
    }
  }
  return (list("N"=N,"N_edges"=N_edges,"node1"=node1,"node2"=node2));
}
nbs <- mungeCARdata4stan(adj, num);

############ get adj matrix 
N <- nbs$N; N_edges <- nbs$N_edges; 
node1 <- nbs$node1; node2 <- nbs$node2;
node_1 <- cbind(node1,node2); node_2 <- cbind(node2,node1)
W <- matrix(0, 83, 83); W[node_1] <- 1; W[node_2] <- 1
W_n <- length(node1); D <- diag(num)
###################Finish getting CAR arguments here######################### 

setwd("..")
save(xx1, xx2, Time, pop_tn, yy, rate, pop_tn_state_plot, teen_d_plot, W, W_n,node1,node2,N_edges, file = "data/collec.RData")

