ESTRELA連載 第78回

今回使うデータ

データの作成と表示

  • パッケージの読み込み
    library(spdep)
    library(gstat)
  • データの読み込み
    # spm.shp <- read.shape("tma_spm.shp")
    spm.shp <- readShapePoints("tma_spm.shp")
    # spm <- spm.shp$att.data
    spm <- cbind(spm.shp$ID, spm.shp$X, spm.shp$Y, spm.shp$SPM07)
    colnames(spm) <- c("ID", "X", "Y", "SPM07")
    spm <- as.data.frame(spm)
    # ward.shp <- read.shape("Ward.shp")
    ward.shp <- readShapePoly("Ward.shp")
    # ward.map <- Map2poly(ward.shp)
    mesh.grid <- read.table("mesh.csv", header=TRUE, sep=",")
    coordinates(mesh.grid) <- c("X", "Y")
    mesh.grid <- as(mesh.grid, "SpatialPixelsDataFrame")

ヴァリオグラム

  • 探索的ヴァリオグラム分析 [#md53f8a0]
  • 等方性モデル
    # 定数項モデル
    spm.var0 <- variogram(SPM07*1000~1, locations=~X+Y, data=spm) 
    plot(spm.var0)
    # 緯度経度によるトレンド
    spm.var1 <- variogram(SPM07*1000~X+Y, locations=~X+Y, data=spm) 
    plot(spm.var1)
    plot(variogram(SPM07*1000~X+Y, locations=~X+Y, data=spm, cloud=TRUE))
  • 異方性モデル
    spm.var2 <- variogram(SPM07*1000~X+Y, locations=~X+Y, data=spm, alpha=0:4*90)
    plot(spm.var2)
  • ヴァリオグラム・モデル
    spm.model1 <- vgm(psill=25, model="Exp", range=28000, nugget=45)
    plot(spm.var1, spm.model1)
    spm.model2 <- vgm(psill=25, model="Sph", range=60000, nugget=45)
    plot(spm.var1, spm.model2)
    spm.model3 <- vgm(psill=25, model="Lin", range=56000, nugget=45)
    plot(spm.var1, spm.model3)
    spm.model4 <- vgm(psill=20, model="Gau", range=35000, nugget=50)
    plot(spm.var1, spm.model4)
    spm.model5 <- vgm(psill=0, model="Nug", nugget=70)
    plot(spm.var1, spm.model5) 
    spm.model6 <- vgm(psill=25, model="Mat", range=30000, nugget=45)
    plot(spm.var1, spm.model6)
  • 異方性モデル
    spm.anis1 <- vgm(psill=25, model="Gau", range=35000, nugget=50, anis=c(0, 0.8))
    plot(spm.var2, spm.anis1)
    # spm.anis2 <- vgm(psill=40, model="Lin", range=60000, nugget=49, anis=c(45, 0.3))
  • ヴァリオグラム・モデルの適合
    spm.fit <- fit.variogram(spm.var1, spm.model1)
    plot(spm.var1, spm.model1)
    # spm.model11 <- vgm(psill=24.55822, model="Exp", range=27627.65, nugget=45.01144)

空間内挿手法

  • 逆距離加重法(IDW)
    spm.idw1 <- idw(SPM07*1000~1, locations=~X+Y, data=spm, mesh.grid, idp=2)
    spplot(spm.idw1["var1.pred"], main = "ordinary kriging predictions")
  • 単純クリギング
    spm.gs <- gstat(id="ID", formula=SPM07*1000~1, locations=~X+Y, data=spm, model=spm.model1, beta=mean(spm$SPM07*1000))
    spm.ps <- predict(spm.gs, mesh.grid)
    spplot(spm.ps[1])
    spplot(spm.ps[2])
  • 通常クリギング
    spm.go <- gstat(id="ID", formula=SPM07*1000~1, locations=~X+Y, data=spm, model=spm.model1)
    spm.po <- predict(spm.go, mesh.grid)
    spplot(spm.po[1])
    spplot(spm.po[2])
  • 普遍クリギング
    spm.gu <- gstat(id="ID", formula=SPM07*1000~X+Y, locations=~X+Y, data=spm, model=spm.model1)
    spm.pu <- predict(spm.gu, mesh.grid)
    spplot(spm.pu[1])
    spplot(spm.pu[2])
  • ブロック・クリギング
    spm.bl <- krige(SPM07*1000~1, locations=~X+Y, data=spm, mesh.grid, model=spm.model1, block=c(50,50))
    spplot(spm.bl[1])
  • シミュレーション:Sequential Gaussian simulation
    spm.sim <- krige(SPM07*1000~1, locations=~X+Y, data=spm, mesh.grid, spm.fit, nsim=6, nmax=100)
    spplot(spm.sim)
  • シミュレーション:Sequential Indicator simulation
    spm.simI <- krige(SPM07*1000~1, locations=~X+Y, data=spm, mesh.grid, spm.fit, indicators = TRUE ,nsim=6, nmax=100)
    spplot(spm.simI)
  • ベイジアン・クリギング

クロス・ヴァリデーション

  • 残差
    mod.set <- sample(1:402,200)
    spm.mod <- spm[mod.set,]
    spm.val <- spm[-mod.set,]
    spm.mod.var <- variogram(SPM07*1000~X+Y, locations=~X+Y, data=spm.mod) 
    plot(spm.mod.var)
    spm.mod.fit <- fit.variogram(spm.mod.var, spm.model1)
    spm.val.pr <- krige(SPM07*1000~X+Y, locations=~X+Y,spm.mod, spm.val, spm.mod.fit)
    spm.res.kr <- spm.val$SPM07*1000 - spm.val.pr$var1.pred
    spm.res.mean <- spm.val$SPM07*1000 - mean(spm.val$SPM07*1000)
    R2 <- 1-sum(spm.res.kr^2)/sum(spm.res.mean^2)
    R2
  • Z値
    spm.cv5 <- krige.cv(SPM07*1000~X+Y, locations=~X+Y, data=spm, spm.fit, nfold=5)

参考

  • クリギング:線形回帰モデル
    # トレンドなし
    spm.krg1 <- krige(SPM07*1000~1, locations=~X+Y, data=spm, mesh.grid, degree=2)
    spplot(spm.krg1["var1.pred"], main = "ordinary kriging predictions")
    # 緯度経度によるトレンド
    spm.krg2 <- krige(SPM07*1000~X+Y, locations=~X+Y, data=spm, mesh.grid)
    spplot(spm.krg2["var1.pred"], main = "ordinary kriging predictions")

トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2010-09-20 (月) 15:10:42 (3343d)