資料處理

library(ggplot2)
library(rgdal)
library(GISTools)
library(sp)
library(aspace)
library(dplyr)
library(spatstat)
library(readxl)
library(spdep)

setwd("D:/spatial/final_project")

pop=readOGR(dsn=".",layer="Popn_TWN2",encoding="utf8")
pop$census=pop$A0A14_CNT+pop$A15A64_CNT+pop$A65UP_CNT
popagg = aggregate(pop$census,by = list(pop$COUNTY),FUN = sum)
names(popagg)[names(popagg)=="Group.1"]="COUNTY"
names(popagg)[names(popagg)=="x"]="census"
hiv=read_excel("hiv.xlsx")
hepa=read_excel("aliver.xlsx")

直方圖：2008年- 2018年的病患人數

#2008-2018 #HIV
hiv_data=data.frame(year=hiv$year,count=as.numeric(as.character(hiv$count)),sex=hiv$sex,nat=hiv$nat) 

agg2008=aggregate(hiv_data$count,by=list(hiv_data$sex,hiv_data$nat,hiv_data$year),FUN=sum)
agg2008=agg2008[21:64,]
#aggdata=data.frame(year=agg2008$Group.1,)

ggplot(agg2008,aes(x=agg2008$Group.3,y=agg2008$x,fill=interaction(Group.2,Group.1,sep = "-",lex.order =TRUE )))+geom_bar(stat="identity")+scale_y_continuous("病患人數")+scale_x_discrete("年份")+labs(title="圖一. 2008～2018HIV患者人數")+scale_fill_manual("國籍與性別",values=c("#F08080","#89CFF0","red1","blue1"))+theme(plot.title=element_text(hjust=0.5))


圖一可以看出在2012年到2017年有增長的趨勢，到了2018年明顯降低。

hepa_data=data.frame(year=hepa$year,count=as.numeric(as.character(hepa$count)),sex=hepa$sex,foreigner=hepa$foreigner)

aggh2008=aggregate(hepa_data$count,by=list(hepa_data$sex,hepa_data$foreigner,hepa_data$year),FUN=sum)
aggh2008=aggh2008[21:64,]
#aggdata=data.frame(year=agg2008$Group.1,)
ggplot(aggh2008,aes(x=aggh2008$Group.3,y=aggh2008$x,fill=interaction(Group.2,Group.1,sep = "-",lex.order =TRUE )))+geom_bar(stat="identity")+scale_y_continuous("病患人數")+scale_x_discrete("年份")+labs(title="圖二. 2008～2018 A型肝炎患者人數")+scale_fill_manual("是否為境外移入 與性別",values=c("#F08080","#89CFF0","red1","blue1"))+theme(plot.title=element_text(hjust=0.5))


從圖二，可以明顯看出在2016年患者人數明顯增加，但到了2017年卻迅速降低。根據資料顯示，這是因為在2015年中，台灣政府開始擴大篩檢範圍，並在家中引入口服液進行A型肝炎檢查。由於從此圖可以看到2016年A型肝炎患者人數為10年內最高，所以接下來的研究會專注於2016年。

面量圖：以病患人數和病患比例分別繪製

hepa=hepa[hepa$year=="2016",] #2016
hepa1=data.frame(count=as.numeric(as.character(hepa$count)),county=as.character(hepa$county),foreigner=hepa$foreigner) 
hepa2016=aggregate(hepa1$count,by=list(hepa1$foreigner,hepa1$county),FUN=sum)

ggplot(hepa2016,aes(x=hepa2016$Group.2,y=hepa2016$x,fill=Group.1))+geom_bar(stat="identity")+scale_y_continuous("病患人數")+scale_x_discrete("縣市")+labs(title="圖三. 2016年各縣市A型肝炎患者人數")+scale_fill_manual("是否為境外移入",values=c("#F08080","#89CFF0"))+theme(plot.title=element_text(hjust=0.5),axis.text.x = element_text(angle=45))


圖四為各縣市的A型肝炎的病例數，可以看到新北市、臺北市和臺中市為前三高的縣市。

hepaagg=aggregate(hepa1$count,by=list(hepa1$county),FUN=sum)
names(hepaagg)[names(hepaagg)=="Group.1"]="COUNTY"
names(hepaagg)[names(hepaagg)=="x"]="patient"
hepaagg1 = full_join(popagg,hepaagg,by="COUNTY")
hepaagg1$propor = hepaagg1$patient/hepaagg1$census*100

popdata=gUnaryUnion(pop,id=pop@data$COUNTY)

row.names(popdata) = as.character(1:length(popdata))
aa=data.frame()
aa=rbind(aa,hepaagg1)
aa=unique(aa$COUNTY)
aa=as.data.frame(aa)
colnames(aa)="COUNTY"

county=SpatialPolygonsDataFrame(popdata,aa)
hepa1=merge(county,hepaagg1,by.x="COUNTY",by.y="COUNTY")
hepapatient=hepa1$patient
shades=auto.shading(hepapatient,n=7,cols=brewer.pal(7,"Greens"))
choropleth(popdata,hepapatient,shades)
choro.legend(379011.7,2531928,shades)
title('圖四. 2016年台灣A型肝炎病患人口分佈圖')
north.arrow(-80122.24,2580664,10000,col="black")
map.scale(-80122.24,2456642,100000,"50km",2,1)

hepapropor=hepa1$propor
shades=auto.shading(hepapropor,n=7,cols=brewer.pal(7,"Greens"))
choropleth(popdata,hepapropor,shades)
choro.legend(379011.7,2531928,shades)
title('圖五. 2016年台灣A型肝炎病患人口比例（%）分佈圖')
north.arrow(-80122.24,2580664,10000,col="black")
map.scale(-80122.24,2456642,100000,"50km",2,1)


從圖四來看，新北市、臺北市、臺中市和桃園市病患人口為最高區間（35-93人）。在圖五，雖然新北市、臺北市和臺中市還是最高層級的，但是桃園市卻相對看起來沒那麼嚴重，高雄市和彰化縣也有了層級上的轉變（變淺）。一個疾病的嚴重度應該取決於人口中患病人數的比例，因為在人口多的地方，患病人數較多為正常現象，因此不可以只是單單看患病人數。

hiv=hiv[hiv$year=="2016",]
hiv1 = data.frame(count=as.numeric(as.character(hiv$count)),county=as.character(hiv$county),nat=hiv$nat)
hiv2016 = aggregate(hiv1$count,by=list(hiv1$nat,hiv1$county),FUN=sum)

ggplot(hiv2016,aes(x=hiv2016$Group.2,y=hiv2016$x,fill=Group.1))+geom_bar(stat="identity")+scale_y_continuous("病患人數")+scale_x_discrete("縣市")+labs(title="圖六. 2016年各縣市HIV患者人數")+scale_fill_manual("國籍",values=c("#F08080","#89CFF0"))+theme(plot.title=element_text(hjust=0.5),axis.text.x = element_text(angle=45))


從圖六，可以發現新北市、臺北市和臺中市依然為患病人數前三高的縣市，因此可以合理推斷A型肝炎和HIV在增長上可能有正向關係。

hivagg=aggregate(hiv$count,by=list(hiv$county),FUN=sum)
names(hivagg)[names(hivagg)=="Group.1"]="COUNTY"
names(hivagg)[names(hivagg)=="x"]="patient"
hivagg1 = full_join(popagg,hivagg,by="COUNTY")
hivagg1$propor = hivagg1$patient/hivagg1$census*100

hiv1=merge(county,hivagg1,by.x="COUNTY",by.y="COUNTY")
hivpatient=hiv1$patient
shades=auto.shading(hivpatient,n=7,cols=brewer.pal(7,"Blues"))
choropleth(popdata,hivpatient,shades)
choro.legend(379011.7,2531928,shades)
title('圖七. 2016年台灣HIV病患人口分佈圖')
north.arrow(-80122.24,2580664,10000,col="black")
map.scale(-80122.24,2456642,100000,"50km",2,1)

hivpropor=hiv1$propor
shades=auto.shading(hivpropor,n=7,cols=brewer.pal(7,"Blues"))
choropleth(popdata,hivpropor,shades)
choro.legend(379011.7,2531928,shades)
title('圖八. 2016年台灣HIV病患人口比例（%）分佈圖')
north.arrow(-80122.24,2580664,10000,col="black")
map.scale(-80122.24,2456642,100000,"50km",2,1)


從兩張圖的比較來看，新北市、臺北市和臺中市還是為最高層級，但是彰化縣和花蓮縣有了明顯的變化。彰化縣變淺了但花蓮縣變深了。花蓮縣只有41個患者人數，相對彰化縣有65個患者人數，從圖上來看，花蓮縣卻看起來比彰化縣嚴重。這是因為花蓮縣人口數明顯較少，所以比例會較高。

皮爾森積差相關分析(Pearson Correlation)：檢測兩種疾病是否存在顯著相關

\({H_0}\):r=0,兩種疾病無顯著相關
\({H_1}\):r \(\neq\) 0,兩種疾病有顯著相關

hiv11 = hiv1$patient
hepa11 = hepa1$patient
plot(hiv11,hepa11,main = "圖九. A型肝炎和HIV的散佈圖",xlab = "HIV", ylab = "A型肝炎",pch = 19)
abline(lm(hepa11~hiv11),col = "red",add = T)
lines(lowess(hiv11,hepa11),col = "blue", add=T)

cor.test(hiv11,hepa11)
## 
##  Pearson's product-moment correlation
## 
## data:  hiv11 and hepa11
## t = 10.809, df = 20, p-value = 8.418e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8230870 0.9683812
## sample estimates:
##       cor 
## 0.9240361


拒絕虛無假設，兩種疾病具有相關性

相關係數為0.92，因此A型肝炎和HIV有很顯著的相關性，這個可以從圖九得到驗證。

繪製LISA 和 Gi* map


由於我們沒有點資料，因此只能使用「面資料」常用區域間彼此自相關的程度來作為地理群聚分析的方法，包括：Anselin’s LISA 及 Getis-Ord’s Gi* (d)。「面資料」為我們將A型肝炎和HIV病例數加總到各縣市，再除以該縣市的人口數作為比例，以各縣市為單位來分析其空間群聚的趨勢。 

coords = coordinates(hepa1)
TW.nb = knn2nb(knearneigh(coords, k=5))
TW.nb.w = nb2listw(TW.nb,zero.policy = T)

LISA = localmoran(hepa1$propor,TW.nb.w,zero.policy = T, alternative = "two.sided" )

propor1 = hepa1$propor
diff = propor1 - mean(propor1) # diff看自己和平均比起來算是H還是L
z = LISA[,4]
quad = c()
quad[diff>0 & z>0] = 1 # H-H
quad[diff<0 & z>0] = 2 # L-L
quad[diff>0 & z<0] = 3 # H-L
quad[diff<0 & z<0] = 4 # L-H
quad[LISA[, 5]>0.1]=5 # 不顯著
colors=c("red", "blue", "lightpink", "skyblue2", rgb(.95, .95, .95))
par(mar=c(0,0,1,0))
plot(hepa1, border="grey", col=colors[quad], main = "圖10. A型肝炎 LISA Map")
legend("bottomright",legend=c("High-High","Low-Low","High-Low","Low-High","No-Sig"),fill=colors,bty="n",cex=0.7,y.intersp=1,x.intersp=1)


圖十是以 Anselin’s LISA 計算A型肝炎局部空間自相關的地理分佈。以最近的前5個縣市設為相鄰定義，圖十的紅色地區表示在0.1的統計顯著水準下，A型肝炎人口比例趨勢顯著地理群聚的區域，可以發現到A型肝炎人口比例在北部地區（新北市、台北市、桃園市）有明顯群聚現象。

TW.nb.in = include.self(TW.nb)
TW.nb.w.in = nb2listw(TW.nb.in)
Gi = localG(propor1,TW.nb.w.in)
LG = as.vector(Gi)

LG1 = 0
for (i in 1:22){LG1[i]<-LG[i]}
hepa1$LG = LG1
lm.palette = colorRampPalette(c("white","orange", "red"), space = "rgb")
spplot(hepa1, zcol="LG", col.regions=lm.palette(20), main="圖十一. A型肝炎 Local Gi*")

# Mapping Hotspot Areas (using Local Gi*, p < 0.05)
quad = c()
quad[LG>=1.645] = 1 # cluster
quad[LG <1.645] = 2 # non-cluster
colors=c("red", "lightgray")
par(mar=c(0,0,1,0))
plot(hepa1, border="grey", col=colors[quad], main = "圖十二. A型肝炎 Cluster Map")
legend("bottomright", c("Cluster","Non-cluster"), fill=colors,bty="n",cex=0.7,y.intersp=1,x.intersp=1)


在圖十一和十二是以Getis-Ord’s Gi\(*\)(d)計算空間自相關的地理分佈。可以看出新竹市、宜蘭縣、桃園市、台北市、新北市和基隆市的A型肝炎人口比例有顯著群聚現象。 Gi\(*\) 在計算鄰居數值時把自己也包括在內。在尋找「熱區」（hotspot）時，把自己的數值也算進去會更接近常理，因為LISA分析的熱區（High-High）是沒有把自己算進去的，所以如果有一個地方自己很高，周圍都很低，就不會出現在熱區裡面。 Gi\(*\)(d) 除了可判定顯著的群聚區域外，更可以呈現其群聚的強度，能夠有效評估疾病風險在地理空間群聚的強度。z值越高，就表示熱點的群聚越緊密，反之z值越小，冷點的聚類就越緊密。 

coords1 = coordinates(hiv1)
TW.nb1 = knn2nb(knearneigh(coords1, k=5))
TW.nb.w1= nb2listw(TW.nb1,zero.policy = T)

LISA1 = localmoran(hiv1$propor,TW.nb.w1,zero.policy = T, alternative = "two.sided" )

propor2 = hiv1$propor
diff1 = propor2 - mean(propor2) # diff看自己和平均比起來算是H還是L
z1 = LISA1[,4]
quad1 = c()
quad1[diff1>0 & z1>0] = 1 # H-H
quad1[diff1<0 & z1>0] = 2 # L-L
quad1[diff1>0 & z1<0] = 3 # H-L
quad1[diff1<0 & z1<0] = 4 # L-H
quad1[LISA1[, 5]>0.1]=5 # 不顯著
colors=c("red", "blue", "lightpink", "skyblue2", rgb(.95, .95, .95))
par(mar=c(0,0,1,0))
plot(hiv1, border="grey", col=colors[quad1], main = "圖十三. HIV LISA Map")
legend("bottomright",legend=c("High-High","Low-Low","High-Low","Low-High","No-Sig"),fill=colors,bty="n",cex=0.7,y.intersp=1,x.intersp=1)

TW.nb.in1 = include.self(TW.nb1)
TW.nb.w.in1 = nb2listw(TW.nb.in1)
Gi1 = localG(propor2,TW.nb.w.in1)
LG2 = as.vector(Gi1)

LG3 = 0
for (i in 1:22){LG3[i]<-LG2[i]}
hiv1$LG = LG3
lm.palette = colorRampPalette(c("white","orange", "red"), space = "rgb")
spplot(hiv1, zcol="LG", col.regions=lm.palette(20), main="圖十四. HIV Local Gi*")

# Mapping Hotspot Areas (using Local Gi*, p < 0.05)
quad1 = c()
quad1[LG2>=1.645] = 1 # cluster
quad1[LG2 <1.645] = 2 # non-cluster
colors=c("red", "lightgray")
par(mar=c(0,0,1,0))
plot(hiv1, border="grey", col=colors[quad1], main = "圖十五. HIV Cluster Map")
legend("bottomright", c("Cluster","Non-cluster"), fill=colors,bty="n",cex=0.7,y.intersp=1,x.intersp=1)


圖十三至十五是以以上的兩種方法，對HIV人口比例進行計算空間自相關的地理分佈。圖十三和十五顯示在桃園市皆有群聚現象。我們認為HIV比較不適合用這種分析群聚的方法來分析，這是因為HIV比較不會透過空間傳染的，但是A型肝炎卻很受空間傳染的影響，如飲用水的來源、人與人之間的接觸等等。