Mid2_score.utf8.md

空間分析　期中考二

助教杜承軒 2020.05.16

圖資讀取

library(rgdal);library(GISTools);library(spatstat);library(splancs);library(magrittr)
setwd("D:/1082SA/Mid2/Data")
case=readOGR("case.shp")
grid=readOGR("grid.shp")
SH=readOGR("SCHOOL.shp",encoding="utf8")
FS=readOGR("FASTFOOD.shp")
TPE=readOGR("TPE.shp")

基礎題

Q1.圖表判讀

[10%]
每一張圖2%*5=10%

b / c / f / g / j

Q2-1. 樣方分析

Q2-1-1
[3%]計算VMR

num=poly.counts(case,grid)
VMR=var(num)/mean(num)
cat("VMR=",VMR,sep='')

VMR=0.5416667

Q2-1-2
[3%]檢定
注意：p-value雙尾(若算成單尾則-1%)

VMR.t=(VMR-1)/sqrt(2/(25-1))
p.two_tail=1-2*abs(0.5-pt(VMR.t,25-1))

p-val=0.1254>0.1，不拒絕虛無假設，病例點呈隨機分布。

Q2-2. 高階最鄰近距離分析

Q2-2-1
[5%]畫圖

cap=ppp(case@coords[,1],case@coords[,2],as.owin(grid))
NNI=colMeans(nndist(cap, k=1:5))
plot(1:5,NNI,type='o',lwd=2,xlab="K",ylab="NND",ylim=c(0,2))

Q2-2-2
[5%]蒙地卡羅檢定
排序後找出前後5個(3%)、畫圖(1%)、結論「不隨機」(1%)

NNI_sim=matrix(nrow=99,ncol=5)
for(i in 1:99) NNI_sim[i,]=colMeans(nndist(rpoint(50, win=as.owin(grid)),k=1:5))
NNI_sim=apply(NNI_sim,2,sort)
plot(1:5,NNI,type='o',lwd=2,xlab="K",ylab="NND",ylim=c(0,2.5))
lines(1:5,NNI_sim[5,],lty=2,type='o')
lines(1:5,NNI_sim[95,],lty=2,type='o')

觀察值絡在信賴包絡之外，拒絕虛無假設，病例點不呈隨機分布。

Q2-3. 高階G函數

Q2-3-1
[4%]畫圖 (若沒有設定k=3則-1%)

ND3=nndist(cap,k=3)
plot(ecdf(ND3),verticals=T,lwd=3,cex=0,col="red",xaxs="i",yaxs ="i",main=expression(G[3](d)),xlab="distance",ylab=expression(G[3](d)))

Q2-3-2
[2%]結論「不隨機」

觀察值絡在信賴包絡之外，拒絕虛無假設，病例點不呈隨機分布。

Q2-4. Ripley’s K函數

[5%]畫圖
扣分：設定nrank=5(-1%)；圖表中L(r)-r，隨機值=0而不是隨機值=r(-1%)；結論「不隨機」(-1%)

L.env=envelope(cap,Lest,nrank=5)
plot(L.env,.-r~r,main="L Function",ylab="L(d)")

觀察值曲線落在信賴包絡之外，拒絕虛無假設，病例點不呈隨機分布。

Q2-5. 綜合分析

[3%]

根據以上檢定方法，評估病例點不為隨機分布。
樣方分析的結果為隨機，建立網格會忽略網格內的點分布，設定網格的大小或位置會影響分析結果。

Q3.程式碼判讀

Q3-1
[4%]

G Function（最鄰近距離的累積頻率分布）。

Q3-2
[6%]
4% + 2%

每個點與最鄰近鄰居的距離。
第一欄是自己和自己距離=0，第二欄才是自己和最近鄰居的距離。

進階題

Q1. 學校肥胖指標差異？

[5%]兩母體平均之t檢定
注意假設和結論要一致（單雙尾）

\(H_0: 公立與私立學校的平均腰圍身高比沒有差異\)
\(H_a: 公立與私立學校的平均腰圍身高比有差異\)
\(\alpha=0.05\)

(T.=t.test(index~Type,SH@data))


    Welch Two Sample t-test

data:  index by Type
t = 5.048, df = 64.98, p-value = 3.849e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.01252157 0.02891541
sample estimates:
mean in group 公立 mean in group 私立 
         0.4349403          0.4142218

p-val<0.05，拒絕虛無假設，公立與私立學校的平均腰圍身高比有差異。

Q2. KDE方法：學校的速食店獨占性

Q2-1
[6%]
※需用splancs的kernel2d，poly引數需從TPE@polygons…擷取
kernel2d設定(2%)、KDE相減(2%)、畫圖(2%)

TPE.bnd=TPE@polygons[[1]]@Polygons[[1]]@coords
SH.KDE=kernel2d(SH@coords,TPE.bnd,1000,100,100,quiet=T)
FS.KDE=kernel2d(FS@coords,TPE.bnd,1000,100,100,quiet=T)
DualKDE=SH.KDE
DualKDE$z=FS.KDE$z-SH.KDE$z
polymap(TPE.bnd,axes=F)
image(DualKDE,add=T,col = brewer.pal(7,"RdYlGn"))
polymap(TPE.bnd,add=T)

Q2-2
[6%]
運算過程(3%)、三間學校的名稱(1%*3=3%)
需先將KDE轉換成SP格式去做套疊，這部分很難，呵呵。
可以用gIntersection的方法應該最直觀，但比較久。
以下用兩個其他方法來作答。

# way 1
  grd=GridTopology(c(DualKDE$x[1],DualKDE$y[1]),c(DualKDE$x[2]-DualKDE$x[1],DualKDE$y[2]-DualKDE$y[1]),c(100,100))
  grd=as.SpatialPolygons.GridTopology(grd,proj4string = TPE@proj4string)
  grd=SpatialPolygonsDataFrame(grd,data=data.frame(z=as.vector(t(apply(DualKDE$z,1,rev))
  )),match.ID=F)
  
  # gIntersection
  # SH.grd=gIntersection(SH,grd,byid=T) 
  # SH.grd.ID=rownames(data.frame(SH.grd))
  # SH.grd.ID=unlist(strsplit(SH.grd.ID," "))
  # grd.ID=SH.grd.ID[seq(2,length(SH.grd.ID),2)]
  # SH$z=grd[grd.ID,]$z
  
  grd.ID=apply(gWithin(SH,grd,byid=T),2,which)
  SH$z=grd[grd.ID,]$z

# way 2
  KDE.sp=SpatialPoints(cbind(rep(DualKDE$x,times=100),rep(DualKDE$y,each=100)),proj4string=FS@proj4string)
  KDE.pts=SpatialPointsDataFrame(KDE.sp,data=data.frame(z=as.vector(DualKDE$z)))
  DIST=gDistance(SH,KDE.pts,byid =T)
  near.id=apply(DIST,2,which.min)
  SH$z=KDE.pts$z[near.id]

High.SH=SH[order(SH$z,decreasing=T)[1:3],]
cat(High.SH$Name%>%as.character())

福星國小 敦化國小 光復國小

附註三間學校的地圖

polymap(TPE.bnd,axes=F)
image(DualKDE,add=T,col = brewer.pal(7,"RdYlGn"))
polymap(TPE.bnd,add=T)
points(High.SH,cex=.8)

Q3. 學校到第K鄰近速食店的距離

Q3-1
[6%]
前十鄰近運算過程(4%)、畫出正確的兩條曲線(2%)

NS=SH[SH$Type=="公立",];PS=SH[SH$Type=="私立",]
NS2FS=rowMeans(apply(gDistance(FS,NS,byid=T),1,sort))
PS2FS=rowMeans(apply(gDistance(FS,PS,byid=T),1,sort))

x=1:10
plot (x,NS2FS[x],type='b',col='blue',xaxt='n',ylim=c(0,3500),xlab="K order",ylab="distance")
lines(x,PS2FS[x],type='b',col='red')
axis(1,x,x)

P.S.也可以用ppp格式+nncross方法，結果圖表一樣。

NP=ppp(NS@coords[,1],NS@coords[,2],as.owin(TPE))
PP=ppp(PS@coords[,1],PS@coords[,2],as.owin(TPE))
FP=ppp(FS@coords[,1],FS@coords[,2],as.owin(TPE))

plot (x,colMeans(nncross(NP,FP,k=x)[,x]),type='b',col='blue',xaxt='n',ylim=c(0,3500),xlab="K order",ylab="distance")
lines(x,colMeans(nncross(PP,FP,k=x)[,x]),type='b',col='red')
axis(1,x,x)

Q3-2
[2%]

學校到鄰近速食店的距離。公立學校離鄰近幾家速食店，平均來說，都比私立學校還近。

Q4. 速食店環形搜尋範圍的肥胖指標

Q4-1
[12%]
運算過程(4%)、速食店的id(2%)、該店指標(2%)、平均指標(4%)。

DIST=gDistance(SH,FS,byid=T)
M=DIST>200&DIST<300
near.index=(M%*%SH$index)/rowSums(M)
max.near=max(near.index,na.rm=T)
max.id=which(near.index==max.near)
mean.index=mean(near.index,na.rm = T)
cat(sprintf("d=300時，鄰近肥胖指標最高的速食店的id為「%d」，指標為「%.4f」。平均鄰近肥胖指標為「%.4f」。",max.id,max.near,mean.index))

d=300時，鄰近肥胖指標最高的速食店的id為「23」，指標為「0.5002」。平均鄰近肥胖指標為「0.4410」。

Q4-2
[6%]
運算過程(4%)、畫出正確的兩條曲線(2%)
P.S.運算過程中應該包含gDistance或gBuffer等相關方法
這小題應該是全部最難的部分

x=seq(100,1000,100)
y.NS=y.PS=c()

DIST.NS=gDistance(NS,FS,byid=T)
DIST.PS=gDistance(PS,FS,byid=T)
for(i in x){
  M.NS=(DIST.NS>=(i-100)&DIST.NS<i)
  M.PS=(DIST.PS>=(i-100)&DIST.PS<i)
  z.NS=M.NS%*%NS$index/rowSums(M.NS)
  z.PS=M.PS%*%PS$index/rowSums(M.PS)
  y.NS=c(y.NS,mean(z.NS,na.rm = T))
  y.PS=c(y.PS,mean(z.PS,na.rm = T))
}

plot(NA,xlim=c(100,1000),ylim=c(0.4,0.53),xaxt='n',xlab="搜尋距離",ylab="平均鄰近肥胖指標")
axis(1,seq(100,1000,100))
lines(x,y.NS,type='o',col='blue')
lines(x,y.PS,type='o',col='red')

Q4-3
[2%]

計算速食店附近學校的肥胖指標，以環形的搜尋區域來觀察，隨著搜尋距離越遠，肥胖指標變化趨勢。能看出隨著距離速食店越遠，平均的肥胖指標有下降再持平的趨勢，且公立的指標皆高於私立。

Q5. 綜合評估

[5%]

速食店、公私立學校與學生肥胖的關連性：
1. 平均而言，公立學校的肥胖指標較私立學校高。
2. 公立學校到速食店距離也比私立學校近。
3. 離速食店越附近的學校，平均的肥胖指標也會較高，且公立學校也高於私立學校。

空間分析 期中考二

助教 杜承軒 2020.05.16

基礎題

進階題

空間分析　期中考二

助教杜承軒 2020.05.16