#obtain acs data
acsdata<-read.csv(file="d:/asc/acs2010-2015.csv",header = TRUE,sep=",")
names(acsdata)
summary(acsdata)
#calculate standard error
acsdata$Poverty_se<-acsdata$Poverty_moe/1.645
acsdata$Edu_lessthan9_se<-acsdata$Edu_lessthan9_moe/1.645
acsdata$Edu_9t12_se<-acsdata$Edu_9t12_moe/1.645
#calculate Education
acsdata$Edu<-acsdata$Edu_9t12+acsdata$Edu_lessthan9
acsdata$Edu_se<-sqrt(acsdata$Edu_lessthan9_se^2+acsdata$Edu_9t12_se^2)
#geting area information
library(tigris)
library(spdep)
counties <- counties(state = 'GA', cb=TRUE)
#create nb list
nb <- poly2nb(counties, queen = FALSE,row.names =counties$GEOID )
attr(nb,"region.id")
results=list()
for (ii in 1:159){
  results[ii]<-nb[which(attr(nb,"region.id")==acsdata$fips[ii])]
}  
dist.n<-matrix(0,159,159)
  for (i in 1:159){
  non.zeroids=unlist(results[i])
  nblist<-attr(nb,"region.id")[non.zeroids]
  non.new<-vector()
  for(ii in 1:length(nblist)){
    non.new[ii]<- which(acsdata$fips==nblist[ii])
  }
  dist.n[i,non.new]=1
  }
#better way
#dist.n<-as.matrix(adj)
adj <- as(dist.n, "dgTMatrix")
#basic setting
betaedu<-0.01
betapov<-0.02
Lambda<-0.05
beta0<-0.1
acsdata$offset<-acsdata$Population*Lambda
acsdata$X<-1:159
#dataset simulation
set.seed(123)

for( ii in 1:200){
  tempdata<-acsdata
  for( i in 1:159){
    tempdata$Edusample[i]<-rnorm(1,acsdata$Edu[i],acsdata$Edu_se[i])
    tempdata$Povertysample[i]<-rnorm(1,acsdata$Poverty[i],acsdata$Poverty_se[i])
    tempdata$y[i]<-rpois(1,acsdata$offset[i]*exp(
      beta0+betaedu*acsdata$Edu[i]+betapov*acsdata$Poverty[i]))
  }
  tempdata$Edusample[tempdata$Edusample<0]<-0
  tempdata$Edusample[tempdata$Edusample>100]<-100
  tempdata$Povertysample[tempdata$Povertysample<0]<-0
  tempdata$Povertysample[tempdata$Povertysample>100]<-100
  assign(paste("data_s", ii, sep=""),tempdata)
}
#fit data to mcmc
library(CARBayes)
formula <- y ~ Edusample+Povertysample+offset(log(offset)) 
for( ii in 1:200){
  tempdata<-eval(parse(text=paste("data_s", ii, sep="")))
  model <- S.CARbym(formula=formula, family="poisson",
                       W=dist.n, burnin=100000, n.sample=200000,data=tempdata)

   assign(paste("mcmcsummary", ii, sep=""),model$summary.results)
   assign(paste("mcmcfit", ii, sep=""),model$fitted.values)
   assign(paste("mcmcdic", ii, sep=""),model$modelfit）
}
#INLA
library(INLA)
for( ii in 1:200){
  tempdata<-eval(parse(text=paste("data_s", ii, sep="")))
  model <- inla(y ~ Edusample+ Povertysample + f(X, model = "bym", graph = adj,
                param=c(1,0.01,1,0.01)) ,control.fixed = list(prec.intercept=0.001),data = tempdata, 
                offset=log(offset),family="poisson",control.predictor = list(compute = TRUE))
  assign(paste("inla", ii, sep=""),model)
}
#simulatio y
sy<-matrix(nrow=200,ncol=159)
for (ii in 1:200){
  tempy<-eval(parse(text=paste("data_s", ii,"$y", sep="")))
  sy[ii,]<-tempy
} 
symean<-colMeans(sy);symean
#rr;
rr<-exp(
  beta0+betaedu*acsdata$Edu+betapov*acsdata$Poverty)
#MCMC results
##beta
mcmcbeta<-data.frame(beta0=double(),
                     betaedu=double(),
                     betapov=double(),
                     tau2=double(),
                     sigma2=double(),
                     stringsAsFactors=FALSE)
for (ii in 1:200){
  tempbeta<-eval(parse(text=paste("mcmcsummary", ii, sep="")))
  mcmcbeta[ii,]<-tempbeta[,1]
} 
mcmcbetamean<-colMeans(mcmcbeta)
mcmcmsebeta<-((beta0-mcmcbetamean[1])^2+(betaedu-mcmcbetamean[2])^2+(betapov-mcmcbetamean[3])^2)/3;mcmcmsebeta


###new
mcmcmse1<-mean((beta0-mcmcbeta[1])^2)
mcmcmse2<-mean((betaedu-mcmcbeta[2])^2)
mcmcmse3<-mean((betapov-mcmcbeta[3])^2)
mcmcmse1
mcmcmse2
mcmcmse2


##relative risk
##Yfitted
mcmcyfit<-matrix(nrow=200,ncol=159)
for (ii in 1:200){
  tempyfit<-eval(parse(text=paste("mcmcfit", ii, sep="")))
  mcmcyfit[ii,]<-tempyfit
} 
mcmcyfitmean<-colMeans(mcmcyfit);mcmcyfitmean
mcmcmsey<-mean((symean-mcmcyfitmean)^2)


####new 
mcmcmsey<-mean(rowMeans((sy-mcmcyfit)^2))
mcmcmsey

#RR
rrmcmc<-cc
mcmcmserr<-mean((rr-rrmcmc)^2)
#new
mcmcrr<-matrix(nrow=200,ncol=159)
for (ii in 1:200){
  mcmcrr[ii,]<-mcmcyfit[1,]/acsdata$offset
} 
mcmcmserr<-mean(rowMeans((rr-mcmcrr)^2))
mcmcmserr

#inla results

##beta
inlabeta<-data.frame(beta0=double(),
                     betaedu=double(),
                     betapov=double(),
                     stringsAsFactors=FALSE)
for (ii in 1:200){
  inlabeta[ii,]<-eval(parse(text=paste("inla", ii,"$summary.fixed[,1]", sep="")))
  
} 
inlabetamean<-colMeans(inlabeta)
inlamsebeta<-((beta0-inlabetamean[1])^2+(betaedu-inlabetamean[2])^2+(betapov-inlabetamean[3])^2)/3;inlamsebeta



##new


inlamse1<-mean((beta0-inlabeta[1])^2)
inlamse2<-mean((betaedu-inlabeta[2])^2)
inlamse3<-mean((betapov-inlabeta[3])^2)
inlamse1
inlamse2
inlamse2
##relative risk
##yfit
inlayfit<-matrix(nrow=200,ncol=159)
for (ii in 1:200){
  tempyfit<-eval(parse(text=paste("inla", ii,"$summary.fitted.values[,1]", sep="")))
  inlayfit[ii,]<-tempyfit
} 
inlayfitmean<-colMeans(inlayfit);inlayfitmean
inlamsey<-mean((symean-inlayfitmean)^2)
####new 
inlamsey<-mean(rowMeans((sy-inlayfit)^2))
inlamsey

##rr
rrinla<-inlayfitmean/acsdata$offset
inlamserr<-mean((rr-rrinla)^2)

#new
inlarr<-matrix(nrow=200,ncol=159)
for (ii in 1:200){
  inlarr[ii,]<-inlayfit[1,]/acsdata$offset
} 
inlamserr<-mean(rowMeans((rr-inlarr)^2))
inlamserr



#y compare
ycompare<-data.frame(x=log(acsdata$Population),
                     ydiffmcmc=symean-mcmcyfitmean,
                     ydiffinla=symean-inlayfitmean,
                     stringsAsFactors=FALSE)
ycompare<-ycompare[order(ycompare$x),]
plot(y=ycompare$ydiffmcmc,x=ycompare$x,type="p",pch=1,col="black",xlab="Log(Population)"
     ,ylab="Difference bettween mean fiited and simulated outcomes"
     ,main="Comparing fit of MCMC and INLA")
lines(y=ycompare$ydiffinla,x=ycompare$x,type="p",pch=4,col="black")
legend(12.5,22,c("MCMC","INLA"),pch=c(1,4),col=c("black","black"))


ycompare1<-data.frame(x=log(acsdata$Population),
                     ydiffmcmc=data_s1$y-mcmcfit1,
                     ydiffinla=data_s1$y-inla1$summary.fitted.values[,1],
                     stringsAsFactors=FALSE)
ycompare1<-ycompare1[order(ycompare1$x),]
plot(y=ycompare1$ydiffmcmc,x=ycompare1$x,type="p",pch=1,col="black",xlab="Log(Population)"
     ,ylab="Difference bettween fiited and simulated outcomes"
     ,main="Comparing fit of MCMC and INLA for simulation dataset 1")
lines(y=ycompare1$ydiffinla,x=ycompare1$x,type="p",pch=4,col="black")
legend(12.5,47,c("MCMC","INLA"),pch=c(1,4),col=c("black","black"))

ycompare2<-data.frame(x=log(acsdata$Population),
                      ydiffmcmc=data_s2$y-mcmcfit2,
                      ydiffinla=data_s2$y-inla2$summary.fitted.values[,1],
                      stringsAsFactors=FALSE)
ycompare2<-ycompare2[order(ycompare2$x),]
plot(y=ycompare2$ydiffmcmc,x=ycompare2$x,type="p",pch=1,col="black",xlab="Log(Population)"
     ,ylab="Difference bettween fiited and simulated outcomes"
     ,main="Comparing fit of MCMC and INLA for simulation dataset 2")
lines(y=ycompare2$ydiffinla,x=ycompare2$x,type="p",pch=4,col="black")
legend(12.5,47,c("MCMC","INLA"),pch=c(1,4),col=c("black","black"))
##map of rr
acsdata$rr<-rr
acsdata$mcmcrr<-rrmcmc
acsdata$inlarr<-rrinla
library(ggplot2)
library(maptools)
library(dplyr)
ggcounty<-fortify(counties, region = "GEOID") 
acsdata$id<-as.character(acsdata$fips)
ggcounty<-left_join(ggcounty, acsdata, by=c("id"))

p1<-ggplot() +
  geom_polygon(data = ggcounty , aes(x=long, y=lat, group = group, fill=Population), color="grey50") +
  scale_fill_gradientn(name="Population",
                       colours = c("black", "white"),
                       values = c(1,0.8, .6, .4, .2, 0))+
  coord_map(xlim = c(-80.5, -86), ylim = c(30,35.5))+
  ggtitle("Population")

p2<-ggplot() +
  geom_polygon(data = ggcounty , aes(x=long, y=lat, group = group, fill=rr), color="grey50") +
  scale_fill_gradientn(name="Relative Risk",
    colours = c("black", "white"),
    values = c(1,0.8, .6, .4, .2, 0))+
  coord_map(xlim = c(-80.5, -86), ylim = c(30,35.5))+
  ggtitle("True Relative Risk")

p3<-ggplot() +
  geom_polygon(data = ggcounty , aes(x=long, y=lat, group = group, fill=mcmcrr), color="grey50") +
  scale_fill_gradientn(name="Relative Risk",
                       colours = c("black", "white"),
                       values = c(1,0.8, .6, .4, .2, 0))+
  coord_map(xlim = c(-80.5, -86), ylim = c(30,35.5))+
  ggtitle("MCMC Relative Risk")

p4<-ggplot() +
  geom_polygon(data = ggcounty , aes(x=long, y=lat, group = group, fill=inlarr), color="grey50") +
  scale_fill_gradientn(name="Relative Risk",
                       colours = c("black", "white"),
                       values = c(1,0.8, .6, .4, .2, 0))+
  coord_map(xlim = c(-80.5, -86), ylim = c(30,35.5))+
  ggtitle("INLA Relative Risk")




multiplot(p1, p2, p3, p4, cols=2)
#locally weighted averages
dist.lwa<-dist.n
for (ii in 1:159){
  dist.lwa[ii,ii]<-1
}  
#edu
edulwa=list()
for (ii in 1:159){
  edulwa[ii]<-sum(dist.lwa[ii,]*acsdata$Edu)/sum(dist.n[ii,])
}  
eduselwa=list()
for (ii in 1:159){
  eduselwa[ii]<-sum(dist.lwa[ii,]*acsdata$Edu_se)/sum(dist.n[ii,])
}  
#pov
povlwa=list()
for (ii in 1:159){
     povlwa[ii]<-sum(dist.lwa[ii,]*acsdata$Poverty)/sum(dist.n[ii,])
}  
povselwa=list()
for (ii in 1:159){
      povselwa[ii]<-sum(dist.lwa[ii,]*acsdata$Poverty_se)/sum(dist.n[ii,])
} 
#combine

acsdata$edulwa<-round(unlist(edulwa),2)
acsdata$eduselwa<-round(unlist(eduselwa),2)
acsdata$povlwa<-round(unlist(povlwa),2)
acsdata$povselwa<-round(unlist(povselwa),2)
#mapping
ggcounty<-fortify(counties, region = "GEOID") 
acsdata$id<-as.character(acsdata$fips)
ggcounty<-left_join(ggcounty, acsdata, by=c("id"))


p1<-ggplot(title="test") +
  geom_polygon(data = ggcounty , aes(x=long, y=lat, group = group, fill=edulwa), color="grey50") +
  scale_fill_gradientn(
    name="%Education",colours = c("black", "white"),
                       values = c(1,0.8, .6, .4, .2, 0))+
  coord_map(xlim = c(-80.5, -86), ylim = c(30,35.5))+
  ggtitle("%Education (Locally Weighted Average)")

p2<-ggplot(title="test") +
  geom_polygon(data = ggcounty , aes(x=long, y=lat, group = group, fill=Edu), color="grey50") +
  scale_fill_gradientn(
    name="%Education",colours = c("black", "white"),
    values = c(1,0.8, .6, .4, .2, 0))+
  coord_map(xlim = c(-80.5, -86), ylim = c(30,35.5))+
  ggtitle("%Education (Raw)")

p3<-ggplot(title="test") +
  geom_polygon(data = ggcounty , aes(x=long, y=lat, group = group, fill=eduselwa), color="grey50") +
  scale_fill_gradientn(name="%SE of Education",colours = c("black", "white"),
                       values = c(1,0.8, .6, .4, .2, 0))+
  coord_map(xlim = c(-80.5, -86), ylim = c(30,35.5))+
  ggtitle("%SE of Education (Locally Weighted Average)")


p4<-ggplot(title="test") +
  geom_polygon(data = ggcounty , aes(x=long, y=lat, group = group, fill=Edu_se), color="grey50") +
  scale_fill_gradientn(name="%SE of Education",colours = c("black", "white"),
                       values = c(1,0.8, .6, .4, .2, 0))+
  coord_map(xlim = c(-80.5, -86), ylim = c(30,35.5))+
  ggtitle("%SE of Education (Raw)")


multiplot(p2, p4, p1, p3, cols=2)

p1<-ggplot(title="test") +
  geom_polygon(data = ggcounty , aes(x=long, y=lat, group = group, fill=povlwa), color="grey50") +
  scale_fill_gradientn(
    name="%Poverty",colours = c("black", "white"),
    values = c(1,0.8, .6, .4, .2, 0))+
  coord_map(xlim = c(-80.5, -86), ylim = c(30,35.5))+
  ggtitle("%Poverty (Locally Weighted Average)")

p2<-ggplot(title="test") +
  geom_polygon(data = ggcounty , aes(x=long, y=lat, group = group, fill=Poverty), color="grey50") +
  scale_fill_gradientn(
    name="%Poverty",colours = c("black", "white"),
    values = c(1,0.8, .6, .4, .2, 0))+
  coord_map(xlim = c(-80.5, -86), ylim = c(30,35.5))+
  ggtitle("%Poverty (Raw)")

p3<-ggplot(title="test") +
  geom_polygon(data = ggcounty , aes(x=long, y=lat, group = group, fill=povselwa), color="grey50") +
  scale_fill_gradientn(name="%SE of Poverty",colours = c("black", "white"),
                       values = c(1,0.8, .6, .4, .2, 0))+
  coord_map(xlim = c(-80.5, -86), ylim = c(30,35.5))+
  ggtitle("%SE of Poverty (Locally Weighted Average)")


p4<-ggplot(title="test") +
  geom_polygon(data = ggcounty , aes(x=long, y=lat, group = group, fill=Poverty_se), color="grey50") +
  scale_fill_gradientn(name="%SE of Poverty",colours = c("black", "white"),
                       values = c(1,0.8, .6, .4, .2, 0))+
  coord_map(xlim = c(-80.5, -86), ylim = c(30,35.5))+
  ggtitle("%SE of Poverty (Raw)")

multiplot(p2, p4, p1, p3, cols=2)
#scrappot 
par(mfrow=c(1,1))
plot(x=acsdata$Edu,y=acsdata$edulwa,xlab="%Education (Raw)",ylab="%Education (Locally Weight Average)")
plot(x=acsdata$Edu_se,y=acsdata$eduselwa,xlab="%SE of Education (Raw)",ylab="%Se of Education (Locally Weight Average)")
plot(x=acsdata$Poverty,y=acsdata$povlwa,xlab="%Poverty (Raw)",ylab="%Poverty (Locally Wegith Average)")
plot(x=acsdata$Poverty_se,y=acsdata$povselwa,xlab="%SE of Poverty (Raw)",ylab="%Se of Poverty (Locally Weight Average)")
#random effect.
inlaran<-matrix(nrow=200,ncol=2)
for (ii in 1:200){
  tempran<-eval(parse(text=paste("inla", ii,"$summary.hyperpar[,1]", sep="")))
  inlaran[ii,]<-tempran
} 
inlaranmean<-1/colMeans(inlaran);inlaranmean

#rrcomparsion
rrcompare<-data.frame(True=rr,
                     mcmc=rrmcmc,
                     inla=rrinla,
                     stringsAsFactors=FALSE)
rrcompare2<-data.frame(y=c(rr,rrmcmc,rrinla),
                      x=c(rep(1,159),rep(2,159),rep(3,159)),
                      stringsAsFactors=FALSE)

plot(x=rrcompare2$x,y=rrcompare2$y,main="Relative Risk Comparsion",ylab="Relative Risk"
     ,xlab="",xaxt = 'n',xlim=c(0.5,3.5))
axis(1, at=1:3, labels=c("Simulation","MCMC","INLA" ))
for (i in 1:159) {
lines(x=rrcompare2$x[c(i,159+i,159*2+i)],y=rrcompare2$y[c(i,159+i,159*2+i)],col="grey")
}



#beta comparsion

mcmcbeta25<-data.frame(beta0=double(),
                     betaedu=double(),
                     betapov=double(),
                     tau2=double(),
                     sigma2=double(),
                     stringsAsFactors=FALSE)
for (ii in 1:200){
  tempbeta<-eval(parse(text=paste("mcmcsummary", ii, sep="")))
  mcmcbeta25[ii,]<-tempbeta[,2]
} 

mcmcbeta975<-data.frame(beta0=double(),
                       betaedu=double(),
                       betapov=double(),
                       tau2=double(),
                       sigma2=double(),
                       stringsAsFactors=FALSE)
for (ii in 1:200){
  tempbeta<-eval(parse(text=paste("mcmcsummary", ii, sep="")))
  mcmcbeta975[ii,]<-tempbeta[,3]
} 

inlabeta25<-data.frame(beta0=double(),
                     betaedu=double(),
                     betapov=double(),
                     stringsAsFactors=FALSE)
for (ii in 1:200){
  inlabeta25[ii,]<-eval(parse(text=paste("inla", ii,"$summary.fixed[,3]", sep="")))
  
} 

inlabeta975<-data.frame(beta0=double(),
                     betaedu=double(),
                     betapov=double(),
                     stringsAsFactors=FALSE)
for (ii in 1:200){
  inlabeta975[ii,]<-eval(parse(text=paste("inla", ii,"$summary.fixed[,5]", sep="")))
  
} 

inlabetase<-data.frame(beta0=double(),
                       betaedu=double(),
                       betapov=double(),
                       stringsAsFactors=FALSE)
for (ii in 1:200){
  inlabetase[ii,]<-eval(parse(text=paste("inla", ii,"$summary.fixed[,2]", sep="")))/sqrt(159)*1.96
  
} 


par(mfrow=c(3,2))

plot(y=1:10,x=mcmcbeta[1:10,1],xlim=c(0.07,0.21),xlab="Intercept",ylab="",yaxt="n",main="Estimated Intercept from MCMC")
abline(v=0.1)
arrows(mcmcbeta25[1:10,1], 1:10, mcmcbeta975[1:10,1],1:10, length=0.05, angle=90, code=3)



plot(y=1:10,x=inlabeta[1:10,1],xlim=c(0.07,0.21),xlab="Intercept",ylab="",yaxt="n",main="Estimated Intercept from INLA")
abline(v=0.1)
arrows(inlabeta25[1:10,1], 1:10, inlabeta975[1:10,1],1:10, length=0.05, angle=90, code=3)



plot(y=1:10,x=mcmcbeta[1:10,2],xlim=c(0.007,0.015),xlab="Education",ylab="",yaxt="n",main="Estimated Education from MCMC")
abline(v=0.01)
arrows(mcmcbeta25[1:10,2], 1:10, mcmcbeta975[1:10,2],1:10, length=0.05, angle=90, code=3)


plot(y=1:10,x=inlabeta[1:10,2],xlim=c(0.007,0.015),xlab="Education",ylab="",yaxt="n",main="Estimated Education from INLA")
abline(v=0.01)
arrows(inlabeta25[1:10,2], 1:10, inlabeta975[1:10,2],1:10, length=0.05, angle=90, code=3)

plot(y=1:10,x=mcmcbeta[1:10,3],xlim=c(0.014,0.021),xlab="Poverty",ylab="",yaxt="n",main="Estimated Poverty from MCMC")
abline(v=0.02)
arrows(mcmcbeta25[1:10,3], 1:10, mcmcbeta975[1:10,3],1:10, length=0.05, angle=90, code=3)


plot(y=1:10,x=inlabeta[1:10,3],xlim=c(0.014,0.021),xlab="Poverty",ylab="",yaxt="n",main="Estimated Poverty from INLA")
abline(v=0.02)
arrows(inlabeta25[1:10,3], 1:10, inlabeta975[1:10,3],1:10, length=0.05, angle=90, code=3)




#mutiple
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the pane
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}


##mcmcconverage
Geweke<-data.frame(beta0=double(),
                     betaedu=double(),
                     betapov=double(),
                     tau2=double(),
                     sigma2=double(),
                     stringsAsFactors=FALSE)
for (ii in 1:200){
  tempbeta<-eval(parse(text=paste("mcmcsummary", ii, sep="")))
  Geweke[ii,]<-tempbeta[,7]
} 
mean(rowMeans(abs(Geweke)<1.96)>0)


#CI%

mean(mcmcbeta25[,1]<=0.1 & mcmcbeta975[,1]>=0.1)
mean(mcmcbeta25[,2]<=0.01 & mcmcbeta975[,2]>=0.01)
mean(mcmcbeta25[,3]<=0.02 & mcmcbeta975[,3]>=0.02)

mean(inlabeta25[,1]<=0.1 & inlabeta975[,1]>=0.1)
mean(inlabeta25[,2]<=0.01 & inlabeta975[,2]>=0.01)
mean(inlabeta25[,3]<=0.02 & inlabeta975[,3]>=0.02)


#moran i
library(spdep)
#edu
col.W <- nb2listw(nb, style="B")
edu <- rep(0,159)

for (i in 1:159) {
 edu[i]<-acsdata$Edu[which(acsdata$fip==as.numeric(attributes(col.W)$region.id[i]))] 
}
str(moran(edu, col.W, length(nb), Szero(col.W)))
#eduse
eduse <- rep(0,159)
for (i in 1:159) {
  eduse[i]<-acsdata$Edu_se[which(acsdata$fip==as.numeric(attributes(col.W)$region.id[i]))] 
}
str(moran(eduse, col.W, length(nb), Szero(col.W)))


#pov
pov <- rep(0,159)

for (i in 1:159) {
  pov[i]<-acsdata$Poverty[which(acsdata$fip==as.numeric(attributes(col.W)$region.id[i]))] 
}
str(moran(pov, col.W, length(nb), Szero(col.W)))
#eduse
povse <- rep(0,159)
for (i in 1:159) {
  povse[i]<-acsdata$Poverty_se[which(acsdata$fip==as.numeric(attributes(col.W)$region.id[i]))] 
}
str(moran(povse, col.W, length(nb), Szero(col.W)))


moran.test(edu,col.W)
moran.test(eduse,col.W)
moran.test(pov,col.W)
moran.test(povse,col.W)

#local
localiedu <- localmoran(edu, col.W)
localieduse <- localmoran(eduse, col.W)
localipov <- localmoran(pov, col.W)
localipovse <- localmoran(povse, col.W)

acsdata$localiedu<-as.numeric(localiedu[,5]<0.05)
acsdata$localieduse<-as.numeric(localieduse[,5]<0.05)

acsdata$localipov<-as.numeric(localipov[,5]<0.05)
acsdata$localipovse<-as.numeric(localipovse[,5]<0.05)
library(spdep)
library(maptools)
library(dplyr)
ggcounty<-fortify(counties, region = "GEOID") 
ggcounty<-left_join(ggcounty, acsdata, by=c("id"))


library(ggplot2)
p2=ggplot(title="test") +
  geom_polygon(data = ggcounty , aes(x=long, y=lat, group = group, fill=factor(localiedu)), color="grey50") +
  scale_fill_manual(name="P-value",
                       values = c(
                         "1"="black",
"0"="white"),label=c("Above or Equal to 0.05","Less than 0.05"))+
  coord_map(xlim = c(-80.5, -86), ylim = c(30,35.5))+
  ggtitle("Local Moran's I for Education")

p1=ggplot(title="test") +
  geom_polygon(data = ggcounty , aes(x=long, y=lat, group = group, fill=factor(localieduse)), color="grey50") +
  scale_fill_manual(name="P-value",
                    values = c(
                      "1"="black",
                      "0"="white"),label=c("Above or Equal to 0.05","Less than 0.05"))+
  coord_map(xlim = c(-80.5, -86), ylim = c(30,35.5))+
  ggtitle("Local Moran's I for SE of Education")

p4=ggplot(title="test") +
  geom_polygon(data = ggcounty , aes(x=long, y=lat, group = group, fill=factor(localipov)), color="grey50") +
  scale_fill_manual(name="P-value",
                    values = c(
                      "1"="black",
                      "0"="white"),label=c("Above or Equal to 0.05","Less than 0.05"))+
  coord_map(xlim = c(-80.5, -86), ylim = c(30,35.5))+
  ggtitle("Local Moran's I for Poverty")

p3=ggplot(title="test") +
  geom_polygon(data = ggcounty , aes(x=long, y=lat, group = group, fill=factor(localipovse)), color="grey50") +
  scale_fill_manual(name="P-value",
                    values = c(
                      "1"="black",
                      "0"="white"),label=c("Above or Equal to 0.05","Less than 0.05"))+
  coord_map(xlim = c(-80.5, -86), ylim = c(30,35.5))+
  ggtitle("Local Moran's I for SE of Poverty")
multiplot(p2, p1, p4, p3, cols=2)