第八章 地理屬性的空間群聚

讀取資料與初步整理

library(sf);library(spdep);library(RColorBrewer);library(units)
setwd("D:/SA_Book/Lab_Data_Ch8")
dengue=st_read(dsn="Dengue_Case.shp")
town=st_read(dsn="KAOH_town.shp")

載入所需之函數

#顏色：紅色為熱區、藍色為冷區，共分成七層
col7=rev(brewer.pal(7,"RdBu")) 

#FDR校正：輸入校正前的Gi*，將輸出FDR校正後的Gi*
FDR.adjust=function(Gi){ 
  Gi=as.numeric(Gi)
  p=mapply(function(x) pnorm(x,lower.tail=x<0),Gi)
  p.fdr=p.adjust(p,"fdr")
  Gi.fdr=qnorm(p.fdr)*-sign(Gi)
  return(Gi.fdr)
}

# 比照ArcGIS PRO中Optimized Hot Spot Analysis功能
# 放入多邊形、數值、鄰近距離閾值三個參數，可得到熱區、冷區顏色
OptimizedHotSpot=function(poly,num,dist){ 
  coord=st_coordinates(st_centroid(poly))
  dnear=dnearneigh(coord, d1=0, d2=dist)
  dnear=include.self(dnear)
  dw=nb2listw(dnear,zero.policy=T,style="W")
  Gi=localG(num,dw)
  Gi.fdr=FDR.adjust(Gi)
  poly$Gi=cut(Gi.fdr,qnorm(c(0,.005,.025,.05,.95,.975,.995,1)))
  poly$Gi_col=col7[poly$Gi]
  return(poly$Gi_col)
}

問題8-1：計算Gi*進行熱區分析

考量到疫情的空間相依特徵，高雄市政府衛生局想進一步找到登革熱病例數量的熱區與冷區，評估各地的嚴重程度。請利用確診病例的點事件資料，以Gi*進行分析。其中，請呈現兩種不同空間尺度的冷熱區分析結果，提供不同地方層級的防疫人員掌握疫情趨勢。

# 六邊形網格
## 建立六邊形網格
hex=st_make_grid(dengue,600,square=F) # 邊長600公尺
hex=st_sf(hex)
## 計算網格內病例數
hex$COUNT=lengths(st_contains(hex,dengue))
hex=hex[hex$COUNT!=0,]
## 熱區冷區著色
hex$col=OptimizedHotSpot(hex,hex$COUNT,6211) #距離參考ArcGIS中參數設定

# 行政區
## 計算行政區內病例數
town$COUNT=lengths(st_contains(town,dengue))
## 熱區冷區著色
town$col=OptimizedHotSpot(town,town$COUNT,15138) #距離參考ArcGIS中參數設定

# 繪圖
par(mar=c(0,0,0,0))
par(mfrow=c(1,2))
plot(st_geometry(hex),col=hex$col)
plot(st_geometry(town),border.col='#DDDDDD',add=T)
plot(st_geometry(town),col=town$col)

問題8-2：進行人口校正的熱區分析

根據問題8-1的分析結果，登革熱感染人數的熱區集中在高雄市南部地區。不過，此處找到的熱區可能是因該地區本身的人口數就較多，故感染人數也較多。若高雄市政府衛生局想要透過冷、熱區來評估地區的危險程度，則須在做局部空間相依性檢測之前，進行人口校正，計算感染人數比例（每多少人中有多少受感者），以調整各地人口數造成的效果。請先針對各區的病例資料進行人口校正，再以Gi*進行分析。

# 計算各行政區的感染人數比例
town$CASE_PRCNT=town$COUNT/town$CENSUS*1000
# 熱區冷區著色
town$col_PRCNT=OptimizedHotSpot(town,town$CASE_PRCNT,15138)

# 繪圖
par(mar=c(0,0,0,0))
plot(st_geometry(town),col=town$col_PRCNT)

問題8-3：最佳化參數設定與統計顯著性校正

在問題8-1與問題8-2中，計算Gi*所使用的Optimized Hot Spot Analysis使用了最佳化的參數設定，包括定義空間關係的帶寬以及使用FDR校正。不過，若視研究主題需要自行調整方法參數，則可使用Hot Spot Analysis (Getis-Ord Gi*) 來計算Gi*。請利用此工具，觀察並比較在檢測局部空間相依性時，重複抽樣會對分析結果造成的影響。

#產生鄰近關係矩陣，此處以距離範圍內來設定鄰近關係，以行政區中心點之間距離來計算
# st_centroid函數：行政區中心點
center=st_centroid(town)
coord=st_coordinates(center)
dnear=dnearneigh(coord, d1=0, d2=14012) #距離參考ArcGIS中參數設定
dnear=include.self(dnear) #Gi*鄰近定義須包含自己
town.dw=nb2listw(dnear,zero.policy=T,style="W") #列標準化

# localG函數：計算Gi*
town$Gi=localG(town$CASE_PRCNT,town.dw)
town$Gi.fdr=FDR.adjust(town$Gi) #FDR校正
# 依照熱區、冷區之90%、95%、99%分成七層，七種顏色
town$Gi_col=col7[cut(town$Gi,qnorm(c(0,.005,.025,.05,.95,.975,.995,1)))] 
town$Gi.fdr_col=col7[cut(town$Gi.fdr,qnorm(c(0,.005,.025,.05,.95,.975,.995,1)))]

# 繪圖
par(mar=c(0,0,0,0))
par(mfrow=c(1,2))
plot(st_geometry(town),col=town$Gi.fdr_col)
plot(st_geometry(town),col=town$Gi_col)

問題8-4：計算Local Moran’s找出熱區與冷區

在分析登革熱感染人數比例局部的空間相依性時，在空間中孤立的冷區代表雖然自己的感染人數比例不高，但四周鄰近地區的疫情嚴重，可能會成為潛在的風險地區。另一方面，如果衛生單位能加強空間中孤立熱區周圍的防疫措施與環境清理，也能夠降低疫情從該區向外擴散的可能性。請透過Local Moran’s I，找到高雄市登革熱疫情孤立的冷區與熱區。

# 產生鄰近關係矩陣，此處以距離矩陣倒數來權重，以行政區中心點之間距離來計算
dist=st_distance(center)
dist=drop_units(dist)
weight=1/dist
weight[dist>14012]=0
diag(weight)=0
town.nbw=mat2listw(weight)

# localmoran_perm函數：計算LISA值(Local Moran's I)
LISA=localmoran_perm(town$CASE_PRCNT,town.nbw,alternative='two.sided')
LISA=as.data.frame(LISA)
LISA$diff=town$CASE_PRCNT-mean(town$CASE_PRCNT)
# 分成H-H/L-L/H-L/L-H/不顯著五組
col=c()
colors=c("red", "blue", "lightpink", "skyblue2", "#F2F2F2")
col[LISA$Z.Ii>0 & LISA$diff>0]=colors[1]
col[LISA$Z.Ii>0 & LISA$diff<0]=colors[2]
col[LISA$Z.Ii<0 & LISA$diff>0]=colors[3]
col[LISA$Z.Ii<0 & LISA$diff<0]=colors[4]
col[LISA$`Pr(z != 0)`>0.1]=colors[5]

# 繪圖
par(mar=c(0,0,0,0))
plot(st_geometry(town),col=col)
legend("bottomright",legend=c("High-High","Low-Low","High-Low","Low-High","Not Significant"), fill=colors,bty="n",cex=0.7,y.intersp=1,x.intersp=1)