空間分析 期末考
助教 杜承軒 2019.06.14
＊初步環境建置與讀取檔案
library(spdep);library(rgdal);library(GISTools);library(magrittr);library(spatstat) setwd("D:/1072SA/Fin") KH=readOGR(dsn =".", layer = "KH", encoding="utf8") TOWN=as.character(KH$TOWN) DG=readOGR(dsn =".", layer = "dengue", encoding="utf8") DG.den=DG$count/poly.areas(DG) TPE=readOGR(dsn =".", layer = "TPE", encoding="utf8")
Q1
請依據以下的鄰近定義，計算a.紅色網格及b.藍色網格分別的統計量，列出計算步驟及結果（列至小數後四位）。
(1)以QUEEN鄰近定義空間鄰近矩陣，手動計算兩個網格的Getis-Ord Local G統計量\(G_i^*\)。[5%]
.=c(2,4,7,9,4,6,3,4,9,4,6,2,6,4,8,1,1,1,7,2,5,10,6,8,6) Gia=sum(2,4,6,3)/sum(.) Gib=sum(4,9,4,6,4,8,1,7,2)/sum(.)(2)以ROOK鄰近並進行列標準化定義空間鄰近矩陣，手動計算兩個網格的Local Moran’s I統計量\(I_i\)。[10%]
.m=mean(.) .v=var(.)*24/25 Iia=(2-.m)*((4-.m)+(6-.m))/.v/2 Iib=(4-.m)*((9-.m)+(6-.m)+(8-.m)+(7-.m))/.v/4(3)透過R建立網格資料，使用相同定義計算上一小題的Local Moran’s I統計量\(I_i\)，並驗證兩者結果是否相同。[10%]
grd=GridTopology(c(0,0),c(1,1),c(5,5)) grd=as.SpatialPolygons.GridTopology(grd) grd=SpatialPolygonsDataFrame(grd,data=data.frame(.),match.ID = F) grd.nb=poly2nb(grd,queen=F) grd.nb.w=nb2listw(grd.nb) grd.LISA=localmoran(grd$.,grd.nb.w) cat(sprintf("(1) Gi* a:% 5.4f , b:% 5.4f \n(2) Ii a:% 5.4f , b:% 5.4f\n(3) Ii a:% 5.4f , b:% 5.4f",Gia,Gib,Iia,Iib,grd.LISA[1,1],grd.LISA[14,1]))## (1) Gi* a: 0.1200 , b: 0.3600 ## (2) Ii a:-0.0000 , b:-0.3634 ## (3) Ii a:-0.0000 , b:-0.3634
Q2
圖資KH.shp為高雄市35個行政區(不含3個原住民區)，其中欄位TOWN及TOWN_ID代表行政區名稱及編號。原本資料內有2018年市長選舉的藍綠比例值數值，透過localmoran函數來檢定哪些行政區的局部空間自相關呈現正相關(顯著水準= 0.05)，但現有的圖資只儲存了計算後原始的Local Moran’s I統計量\(I_i\)和其期望值\(E(I_i)\)、變異數\(Var(I_i)\)。(1)計算p-value，並列出局部空間自相關呈現正相關(α= 0.05)的行政區名稱。[5%]
zi=(KH$Ii-KH$E_Ii)/sqrt(KH$Var_Ii) pi=pnorm(zi,lower.tail = F) TOWN[pi<0.05]## [1] "左營區" "美濃區" "甲仙區" "杉林區"(2)對所有行政區透過手動計算方式來進行FDR校正，由小到大依序列出校正後p-value前三低的行政區名稱。[10%]
p.fdr=pi*length(pi)/rank(pi) p.fdr[p.fdr>1]=1 TOWN[order(p.fdr)[1:3]]## [1] "杉林區" "甲仙區" "左營區"(3)透過p.adjustSP函數，考慮QUEEN的鄰近關係來校正FDR，並列出局部空間自相關呈現正相關(α= 0.05)的行政區名稱。[10%]
KH.nb=poly2nb(KH) p.fdr.sp=p.adjustSP(pi,KH.nb,method='fdr') TOWN[p.fdr.sp<0.05]## [1] "左營區" "甲仙區" "杉林區"(4)上述(2)、(3)小題都是在進行FDR校正，請以文字說明，兩者概念與實作的差異為何？並比較那一個方法得出的結果較保守。[5%]
cat('前者以所有多邊形個數計算，後者考慮鄰居的個數來進行校正。前者較保守。')## 前者以所有多邊形個數計算，後者考慮鄰居的個數來進行校正。前者較保守。
Q3
圖資dengue.shp為某年南高屏三縣市登革熱爆發的病例，其中欄位count為病例數。請以各鄉鎮病例密度（即每單位面積有多少病例）作為參考數值，以距離作為鄰近基準，來量測空間自相關。(1)評估各鄉鎮中心之間的距離，決定鄰近距離d，並說明選擇的原因。[5%]
coord=coordinates(DG) x=coord[,1];y=coord[,2] cen=ppp(x,y,range(x),range(y)) plot(c(0,20000),c(0,1),'n',xlab='距離',ylab='F 函數',xaxs = "i", yaxs ="i") for(i in 1:8) lines(ecdf(nndist(cen,k=i)),col=i,cex=0,verticals=T)
d=10000(2)以上述(1)的距離，以一階鄰近(距離0~d)為定義，繪製出LISA地圖(以顯著水準=0.1，呈現出HH、LL、HL、LH)。[5%]
LISA=function(diff,lisa,p){ z=lisa[,4] quad=c() quad[diff>0 & z>0]=1;quad[diff<0 & z>0]=2;quad[diff>0 & z<0]=3;quad[diff<0 & z<0]=4;quad[lisa[, 5]>p]=5 return(quad) } LISA.legend=function(){legend("bottomright",legend=c("HH","LL","HL","LH","Not-Sig"), fill=colors,bty="n",cex=0.7,y.intersp=1,x.intersp=1)} colors=c("red", "blue", "lightpink", "skyblue2","#EEEEEE") TW.d=dnearneigh(coord,d1=0, d2=d) TW.dw=nb2listw(TW.d,zero.policy=T) DG.d1=localmoran(DG.den,TW.dw,zero.policy=T,alternative="two.sided") quad=LISA(DG.den-mean(DG.den),DG.d1,0.1) plot(DG, border="grey", col=colors[quad]) LISA.legend()
(3)以第二階鄰近(距離d~2d)為定義，繪製出LISA地圖(以顯著水準=0.1，呈現出HH、LL、HL、LH)。[10%]
TW.d2=dnearneigh(coord,d1=d, d2=d*2) TW.dw2=nb2listw(TW.d2,zero.policy=T) DG.d2=localmoran(DG.den,TW.dw2,zero.policy=T,alternative="two.sided") quad=LISA(DG.den-mean(DG.den),DG.d2,0.1) plot(DG, border="grey", col=colors[quad]) LISA.legend()
(4)說明(2)、(3)地圖代表的地理意涵。[5%]
cat('第一階鄰近捕捉熱區，第二階鄰近捕捉登革熱傳播隨距離遞減的效應')## 第一階鄰近捕捉熱區，第二階鄰近捕捉登革熱傳播隨距離遞減的效應
Q4
圖資TPE.shp為台北市次分區，其中欄位DEN代表個次分區中的人口密度。以QUEEN鄰近並進行列標準化定義鄰近關係，來計算空間自相關的數值。過去課程實習在檢定的方法，都是以moran.mc或localmoran函數中已經執行的結果，如下圖。
希望能用手動模擬的方式，也能得到類似的結果，因此，請透過排列檢定(permutation test)，模擬99次的隨機情形，比較觀察值與隨機值。
(1)繪製出Global Moran’s I的直方圖。[10%]
(2) 透過與隨機值的排序比較，計算各個區局部空間自相關呈現正相關的p-value。並以data.frame的方式，呈現行政區名、localmoran函數算出來的p-value、模擬計算p-value，如下圖。[10%]
(注意：模擬次數加上觀察值為100次，故模擬計算p-value一定是二位以下的小數；而下圖localmoran函數算出來的p-value有經過四捨五入)## TPE p_org p_sim ## 0 三民地區 0.60 0.64 ## 1 東社地區 0.63 0.69 ## 2 中崙地區 0.04 0.02 ## 3 本鎮地區 0.04 0.01 ## 4 三張犁地區 0.08 0.05TW.c=poly2nb(TPE) TW.cw=nb2listw(TW.c) den=TPE$DEN mc=moran.mc(den,listw=TW.cw,nsim=99) MoranI=mc$statistic LISA=localmoran(den,listw=TW.cw) Ii=LISA[,1] for(i in 1:99) { s=sample(den) MoranI=c(MoranI,moran.test(s,listw=TW.cw)$estimate[1]) Ii=cbind(Ii,localmoran(s,listw=TW.cw)[,1]) } hist(MoranI, freq=TRUE, breaks=10, xlab="Simulated Moran's I",main="Monte-Carlo simulation") abline(v=MoranI[1], col="red")
p.sim=apply(Ii,1,function(x) as.numeric((101-rank(x))[1])/100) data.frame(TPE=TPE$SEC_NAME,p_org=round(LISA[,5],2),p.sim)## TPE p_org p.sim ## 0 三民地區 0.60 0.68 ## 1 東社地區 0.63 0.77 ## 2 中崙地區 0.04 0.09 ## 3 本鎮地區 0.04 0.07 ## 4 三張犁地區 0.08 0.07 ## 5 五分埔地區 0.10 0.08 ## 6 福德地區 0.47 0.42 ## 7 吳興地區 0.26 0.30 ## 8 六張犁地區 0.45 0.39 ## 9 和平地區 0.20 0.12 ## 10 學府地區 0.63 0.69 ## 11 瑞安地區 0.00 0.01 ## 12 新生地區 0.19 0.11 ## 13 敦南地區 0.00 0.01 ## 14 安和地區 0.00 0.01 ## 15 臥龍地區 0.42 0.28 ## 16 林森地區 0.31 0.23 ## 17 圓山地區 0.25 0.12 ## 18 新庄地區 0.41 0.37 ## 19 大直地區 0.23 0.18 ## 20 下埤頭地區 0.52 0.58 ## 21 長春地區 0.26 0.16 ## 22 朱厝崙地區 0.93 0.88 ## 23 古亭地區 0.26 0.17 ## 24 公館地區 0.69 0.76 ## 25 南門地區 0.00 0.01 ## 26 崁頂地區 0.00 0.02 ## 27 城內地區 0.38 0.30 ## 28 東門地區 0.37 0.20 ## 29 大龍地區 0.44 0.38 ## 30 蘭州地區 0.48 0.48 ## 31 建成地區 0.16 0.10 ## 32 延平地區 0.40 0.26 ## 33 西門地區 0.46 0.43 ## 34 龍山地區 0.52 0.62 ## 35 大理地區 0.02 0.07 ## 36 西園地區 0.34 0.22 ## 37 東園地區 0.80 0.84 ## 38 青年地區 0.06 0.08 ## 39 景美地區 0.41 0.31 ## 40 興隆地區 0.47 0.40 ## 41 木柵地區 0.32 0.25 ## 42 二格山地區 0.06 0.08 ## 43 萬芳地區 0.07 0.08 ## 44 三重埔地區 0.03 0.08 ## 45 新庄仔地區 0.07 0.09 ## 46 後山坡地區 0.47 0.50 ## 47 中研地區 0.00 0.03 ## 48 東湖地區 0.04 0.06 ## 49 西湖地區 0.06 0.08 ## 50 紫陽地區 0.94 0.95 ## 51 金龍地區 0.03 0.03 ## 52 灣仔地區 0.20 0.15 ## 53 洲尾地區 0.49 0.50 ## 54 街上地區 0.09 0.06 ## 55 後港地區 0.52 0.64 ## 56 社子地區 0.27 0.16 ## 57 芝山岩地區 0.15 0.12 ## 58 蘭雅地區 0.85 0.97 ## 59 天母地區 0.18 0.12 ## 60 陽明山地區 0.00 0.01 ## 61 石牌地區 0.09 0.07 ## 62 唭哩岸地區 0.45 0.48 ## 63 關渡地區 0.00 0.02 ## 64 舊北投地區 0.15 0.11 ## 65 新北投地區 0.12 0.10 ## 66 陽明山地區 0.00 0.01 ## 67 大屯地區 0.00 0.02
