環境ガバナンスのデータサイエンス/空間モデリング特論　R演習7

演習の概要

• カーネル密度関数、空間内挿

注意事項

• 慶應義塾大学SFCで開講している「環境ガバナンスのデータサイエンス」「空間モデリング特論」の授業履修者向けの演習用ページです。
• 必ずしも全てのバージョンのRやOSで動作確認を行っていません。この演習用ページを作成している段階では、R3.5.3を使っています。
• Rの更新などにより、Rコードなどを予告無しに変更する場合があります。

データ分析の準備

データのダウンロード

• 演習で用いるデータはここからダウンロードしてください。
• ここでは./直下にgisdataフォルダを作成し、setwd("gisdata")とディレクトリを指定しています

パッケージのインストール

• install.packages()でパッケージをインストールし、library()でパッケージを呼び出す
• ここでは以下のパッケージを用います

install.packages("spdep")
install.packages("spatstat")
install.packages("splancs")
install.packages("spatialkernel")
install.packages("spatial")
install.packages("gstat", dependencies = T)
install.packages("sp")

library(spdep)
library(spatstat)
library(splancs)
library(spatialkernel)
library(spatial)
library(gstat)
library(sp)

以下の方法でパイプ %>% の優先順位を高める

needs::prioritize(magrittr)

空間内挿

カーネル密度関数

まず、カーネル密度関数の基本となる、ガウス関数、イパクニコフ関数、四次関数を図示してみよう。 #### ガウス関数

plot(density(0, kernel="gaussian", bw=1), 
     main="", xlab="", ylab="",lwd=2, cex.axis=1)

イパクニコフ関数

plot(density(0, kernel="epanechnikov", bw=1), 
     main="", xlab="", ylab="",lwd=2, cex.axis=1)

四次関数

plot(density(0, kernel="biweight", bw=1), 
     main="", xlab="", ylab="",lwd=2, cex.axis=1)

バンド幅とカーネル密度関数（合成関数）との関係

バンド幅を変えた場合にカーネル密度関数の合成関数がどのように変化するかを見てみよう。

ここでは、ガウス関数を用い、バンド幅=1とした場合の例を示す。

# (a) バンド幅=1
ss <-1
# ss <-0.7    (b) バンド幅=0.7
# ss <-2      (c) バンド幅=2
gau.all <- function(x){
1/sqrt(2*pi*ss)*exp(-x^2/(2*ss^2)) + 
1/sqrt(2*pi*ss)*exp(-(x+4)^2/(2*ss^2)) +
1/sqrt(2*pi*ss)*exp(-(x-2)^2/(2*ss^2)) +
1/sqrt(2*pi*ss)*exp(-(x-5)^2/(2*ss^2)) +
1/sqrt(2*pi*ss)*exp(-(x+1.5)^2/(2*ss^2))}
curve(gau.all, xlim=c(-10, 10), ylim=c(0, 0.8), main="", xlab="", ylab="",
lwd=8, cex.axis=1.8)
#
gau1 <-  function(x){1/sqrt(2*pi*ss)*exp(-x^2/(2*ss^2))}
gau2 <-  function(x){1/sqrt(2*pi*ss)*exp(-(x+4)^2/(2*ss^2))}
gau3 <-  function(x){1/sqrt(2*pi*ss)*exp(-(x-2)^2/(2*ss^2))}
gau4 <-  function(x){1/sqrt(2*pi*ss)*exp(-(x-5)^2/(2*ss^2))}
gau5 <-  function(x){1/sqrt(2*pi*ss)*exp(-(x+1.5)^2/(2*ss^2))}
#
curve(gau1, lwd=4, add=TRUE)

curve(gau2, lwd=4, add=TRUE)
curve(gau3, lwd=4, add=TRUE)

curve(gau4, lwd=4, add=TRUE)
curve(gau5, lwd=4, add=TRUE)

２次元平面空間上での密度関数

次に、２次元平面空間上でバンド幅を変えた場合の影響を見てみよう。

ここでは、ガウス関数を用い、バンド幅＝１年経例を示す。

x <- runif(50)*10
y <- runif(50)*10
xy <- cbind(x, y)
poly0 <- cbind(c(0,10,10,0), c(0,0,10,10))
# バンド幅＝１
image(splancs::kernel2d(xy, poly0, h0=1, nx=100, ny=100),
      col=gray((0:20)/20),cex.axis=1.8)

Xrange is  0 10 
Yrange is  0 10 
Doing quartic kernel

# バンド幅＝0.5
# image(kernel2d(xy, poly0, h0=0.5, nx=100, ny=100), col=gray((0:20)/20),cex.axis=1.5)
# バンド幅＝0.7
# image(kernel2d(xy, poly0, h0=0.7, nx=100, ny=100), col=gray((0:20)/20),cex.axis=1.5)

最小二乗法によるバンド幅の決定

splancsパッケージのmse2d()関数を用いて、最小二乗法によりバンド幅を決定することができる。

poly0 <- cbind(c(0,1,1,0), c(0,0,1,1))
X <- cbind(runif(50), runif(50))
Mse2d <- splancs::mse2d(splancs::as.points(X), poly0, nsmse=50, range=1)
plot(Mse2d$h[5:50],Mse2d$mse[5:50], type="l", ylab="MSE", xlab="バンド幅",
     lwd=1, cex.axis=1, cex.lab=1)
points(Mse2d$h[which.min(Mse2d$mse)], Mse2d$mse[which.min(Mse2d$mse)],
pch=1, cex=3, lwd=2)

交差検証対数尤度関数によるバンド幅の決定

spatialkernelパッケージのcvloglk() 関数を用いて、交差検証対数尤度関数によるバンド幅を決定することができる。

x <- runif(300)
y <- runif(300)
mks <- sample(c("a","b", "c"), 300, replace=TRUE)
pts <- cbind(x, y)
h <- seq(0.01, 1, by=0.01)
cv <- spatialkernel::cvloglk(pts, mks, h=h)$cv
# 図9.6
plot(h, cv, type="l", ylab="交差検証対数尤度", xlab="バンド幅",
     lwd=6, cex.axis=1.3, cex.lab=1.2)
points(h[which.max(cv)], cv[which.max(cv)], pch=1, cex=3, lwd=2)

データの読み込み

spm.shp  <- sf::st_read('tma_spm.shp')
spm <- ppp(spm.shp$X, spm.shp$Y, c(-70000, 90000), c(-110000, 30000),
           marks=spm.shp$SPM07*1000)
ward.shp <- sf::st_read('Ward.shp')
mesh.grid <- readr::read_csv('mesh.csv')
coordinates(mesh.grid) <- c("X", "Y")
mesh.grid <- as(mesh.grid, "SpatialPixelsDataFrame")

逆距離加重法

spm.idw1 <- idw(X=spm, power=2, at="pixels", se=TRUE)
plot(spm.idw1$estimate)

plot(spm.idw1$SE)

クリギング

データの準備

spm.v <- cbind(spm.shp$ID, spm.shp$X, spm.shp$Y, spm.shp$SPM07)
colnames(spm.v) <- c("ID", "X", "Y", "SPM07")
spm.v <- as.data.frame(spm.v)
sp::coordinates(spm.v) <- c("X","Y")

ヴァリオグラム・モデル

ヴァリオグラム雲と標本ヴァリオグラム

まずは定数項を考慮したモデルを推定し、ヴァリオグラム雲と標本ヴァリオグラムを可視化してみよう。

spm.var1 <- variogram(SPM07*1000~X+Y, data=spm.v) 

variogram()関数のcloud=TRUEとすると、ヴァリオグラム雲を作成できる。

var.cld <- variogram(SPM07*1000~X+Y, data=spm.v, cloud=TRUE)
plot(var.cld$dist, var.cld$gamma, pch=19, cex=1, 
  xlab="distance", ylab="gamma", cex.axis=1, cex.lab=1)

また標本ヴァリオグラムは以下のように可視化できる。

plot(spm.var1$dist, spm.var1$gamma, pch=1, lwd=1, cex=1, ylim=c(0, 80),
  xlab="distance", ylab="gamma", cex.axis=1, cex.lab=1)

指数モデル

plot(variogramLine(vgm(psill=25, model="Exp", range=28000, nugget=45),
70000, 100), type="l", cex=1.5, lwd=4, ylab="semivariance", xlab="distance",
cex.axis=1.2, cex.lab=1.2, ylim=c(40,80))
points(spm.var1$dist, spm.var1$gamma, cex=1.5, lwd=2) 

#spm.model1 <- vgm(psill=25, model="Exp", range=28000, nugget=45)
#plot(spm.var1, spm.model1, cex=1.5, lwd=4)

球形モデル

plot(variogramLine(vgm(psill=25, model="Sph", range=60000, nugget=45),
70000, 100), type="l", cex=1.5, lwd=4, ylab="semivariance", xlab="distance",
cex.axis=1.2, cex.lab=1.2, ylim=c(40,80))
points(spm.var1$dist, spm.var1$gamma, cex=1.5, lwd=2) 

# spm.model2 <- vgm(psill=25, model="Sph", range=60000, nugget=45)
# plot(spm.var1, spm.model2, cex=1.5, lwd=4)

線形モデル

plot(variogramLine(vgm(psill=25, model="Lin", range=56000, nugget=45),
70000, 100), type="l", cex=1.5, lwd=4, ylab="semivariance", xlab="distance",
cex.axis=1.2, cex.lab=1.2, ylim=c(40,80))
points(spm.var1$dist, spm.var1$gamma, cex=1.5, lwd=2) 

# spm.model3 <- vgm(psill=25, model="Lin", range=56000, nugget=45)
# plot(spm.var1, spm.model3, cex=1.5, lwd=4)

ガウスモデル

plot(variogramLine(vgm(psill=20, model="Gau", range=35000, nugget=50),
70000, 100), type="l", cex=1.5, lwd=4, ylab="semivariance", xlab="distance",
cex.axis=1.2, cex.lab=1.2, ylim=c(40,80))
points(spm.var1$dist, spm.var1$gamma, cex=1.5, lwd=2) 

# spm.model4 <- vgm(psill=20, model="Gau", range=35000, nugget=50)
# plot(spm.var1, spm.model4, cex=1.5, lwd=4)

ナゲット効果モデル

plot(variogramLine(vgm(psill=0, model="Nug", range=0, nugget=70),
70000, 100), type="l", cex=1.5, lwd=4, ylab="semivariance", xlab="distance",
cex.axis=1.2, cex.lab=1.2, ylim=c(40,80))
points(spm.var1$dist, spm.var1$gamma, cex=1.5, lwd=2) 

# spm.model5 <- vgm(psill=0, model="Nug", nugget=70)
# plot(spm.var1, spm.model5, cex=1.5, lwd=4) 

Maternモデル

plot(variogramLine(vgm(psill=25, model="Mat", range=30000, nugget=45),
70000, 100), type="l", cex=1.5, lwd=4, ylab="semivariance", xlab="distance",
cex.axis=1.2, cex.lab=1.2, ylim=c(40,80))
points(spm.var1$dist, spm.var1$gamma, cex=1.5, lwd=2) 

# spm.model6 <- vgm(psill=25, model="Mat", range=30000, nugget=45)
# plot(spm.var1, spm.model6, cex=1.5, lwd=4)

推定方法によるヴァリオグラム・モデルの違い

ヴァリオグラム・モデルは推定方法によって推定結果が異なる場合がある。ここでは、球形モデルを例に、重み付き最小二乗法、通常最小二乗法、制限付き最尤法の３つの推定方法による推定結果の違いを示そう。

spm.model2 <- vgm(psill=25, model="Sph", range=60000, nugget=45)
spm.fit <- fit.variogram(spm.var1, spm.model2)
# 球形モデル：WLS
fit.variogram(spm.var1, spm.model2, fit.method=7)

# 球形モデル：OLS
fit.variogram(spm.var1, spm.model2, fit.method=6)

# 制限付き最尤法
fit.variogram.reml(SPM07*1000~X+Y,data=spm.v, model=vgm(25, "Sph", 60000, 45))

# WLS
plot(variogramLine(vgm(psill=25, model="Sph", range=60000, nugget=45),
70000, 100), type="l", cex=1.5, lwd=4, ylab="semivariance", xlab="distance",
cex.axis=1.2, cex.lab=1.2, ylim=c(40,80))
# OLS
lines(variogramLine(vgm(psill=25.15494, model="Sph", range=93784.44, nugget=48.34058), 100000, 100), cex=1.5, lwd=4,lty=2)

# REML
lines(variogramLine(vgm(psill=28.94858, model="Sph", range=60000, nugget=51.51024), 100000, 100), cex=1.5, lwd=4,lty=3)
# バリオグラム
points(spm.var1$dist, spm.var1$gamma, cex=1.5, lwd=2) 

# 凡例
legend("bottomright", legend=c("WLS", "OLS", "REML"), lty=c(1,2,3), lwd=c(4,4,4), cex=1.5)

異方性モデリング

variogram()関数のalphaを指定すると、異方性(anisotropy)を考慮したヴァリオグラム・モデルを推定できる。

spm.var2 <- variogram(SPM07*1000~X+Y, data=spm.v, alpha=0:3*90)
plot(spm.var2)

また、anisを指定して異方性モデルを推定することもできる。

spm.anis1 <- vgm(psill=25, model="Gau", range=35000, nugget=50, anis=c(0, 0.8))
plot(spm.var2, spm.anis1)

クリギングによる内挿補間

ヴァリオグラム・モデルに指数モデルを採用し、クリギングに空間内挿補間を適用しよう。

ここでは、単純クリギング、通常クリギング、普遍クリギング、ブロック・クリギングの適用例を示す。 ##### 単純クリギング

#spm.model1 <- vgm(psill=25, model="Exp", range=28000, nugget=45)
spm.krige1 <- krige(SPM07*1000~1, spm.v, mesh.grid, model=spm.model1,
                   beta=mean(spm.v$SPM07*1000))

[using simple kriging]

spplot(spm.krige1["var1.pred"], main="Simple Kriging Prediction")

spplot(spm.krige1["var1.var"], main="Simple Kriging Variance")

通常クリギング

spm.krige2 <- krige(SPM07*1000~1, spm.v, mesh.grid, model=spm.model1)

[using ordinary kriging]

spplot(spm.krige2["var1.pred"], main="Ordinary Kriging Prediction")

spplot(spm.krige2["var1.var"], main="Ordinary Kriging Variance")

普遍クリギング

spm.krige3 <- krige(SPM07*1000~X+Y, spm.v, mesh.grid, model=spm.model1, 
                    block=c(250,250))

[using universal kriging]

spplot(spm.krige3["var1.pred"], main="Universal Kriging Prediction")

spplot(spm.krige3["var1.var"], main="Universal Kriging Variance")