sdmTMBによる空間変量効果を組み込んだモデル

R
TMB
作者

伊東宏樹

公開

2024年9月20日

TMBを使って空間分布モデルを扱うRパッケージにsdmTMBがあります。今回はsdmTMBをつかって、空間変量効果を組み込んだモデルのあてはめをおこなってみました。

準備

sdmTMBパッケージと、グラフ表示のためggplot2パッケージを読み込みます。

library(sdmTMB)
library(readr)
library(ggplot2)

データ

今回も、前回に続いて、岩波データサイエンスVol.1で使用した、Y方向10×X方向20のグリッド中のアラカシの株数データを使用します。

qgl <- read_csv(url("https://raw.githubusercontent.com/iwanami-datascience/vol1/master/ito/Qglauca.csv"))  |>
  dplyr::arrange(Y, X)
## Rows: 200 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (3): X, Y, N
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(qgl)
## # A tibble: 6 × 3
##       X     Y     N
##   <dbl> <dbl> <dbl>
## 1     1     1     0
## 2     2     1     0
## 3     3     1     0
## 4     4     1     0
## 5     5     1     0
## 6     6     1     0

図示

いつものように、まずデータを図示します。

ggplot(qgl, aes(x = X, y = Y)) +
  geom_tile(aes(fill = N)) +
  scale_y_reverse(breaks = c(0, 5, 10)) +
  scale_fill_gradient(name = "株数", low = "grey90", high = "grey10",
                      limits = c(0, 6)) +
  coord_fixed(ratio = 1) +
  theme_bw(base_family = "IPAexGothic")

メッシュの作成

sdmTMBでは、空間データはメッシュで与えることになっています。そのためのメッシュをmake_mesh関数で作成します。ここではtype = "cutoff"としていますが、typeには"kmeans"もえらべます。

まずは、メッシュの辺の最小の長さが2となるようにcutoff = 2としています。plotで作成されたメッシュを確認できます。

mesh1 <- make_mesh(qgl, xy_cols = c("X", "Y"), type = "cutoff", cutoff = 2)
plot(mesh1)

あてはめ

sdmTMB関数でモデルをあてはめます。ここではY座標の値を説明変数として、空間変量効果をふくめ、誤差構造をポアソン分布、リンク関数をlogとするモデルとしています。

fit1 <- sdmTMB(N ~ Y, data = qgl, family = poisson(link = "log"),
               mesh = mesh1, spatial = "on")
print(fit1)
## Spatial model fit by ML ['sdmTMB']
## Formula: N ~ Y
## Mesh: mesh1 (isotropic covariance)
## Data: qgl
## Family: poisson(link = 'log')
##  
##             coef.est coef.se
## (Intercept)    -1.63    0.65
## Y               0.21    0.09
## 
## Matérn range: 7.98
## Spatial SD: 0.66
## ML criterion at convergence: 259.406
## 
## See ?tidy.sdmTMB to extract these values as a data frame.

Yの係数は0.21と推定されました。

予測値

予測値を図示します。

p1 <- predict(fit1, newdata = dplyr::select(qgl, X, Y))
qgl$est1 <- exp(p1$est)
ggplot(qgl, aes(x = X, y = Y)) +
  geom_tile(aes(fill = est1)) +
  scale_y_reverse(breaks = c(0, 5, 10)) +
  scale_fill_gradient(name = "株数", low = "grey90", high = "grey10",
                      limits = c(0, 6)) +
  coord_fixed(ratio = 1) +
  theme_bw(base_family = "IPAexGothic")

わりとなめらかに変化する予測値が得られました。

メッシュを変えた場合

次にメッシュの大きさを変えて、予測値がどのように変化するか見てみます。

cutoff = 0.9として、メッシュを作成します。

mesh2 <- make_mesh(qgl, xy_cols = c("X", "Y"), type = "cutoff", cutoff = 0.9)
plot(mesh2)

メッシュが細かくなりました。

あてはめ

モデルをあてはめます。

fit2 <- sdmTMB(N ~ Y, data = qgl,family = poisson(link = "log"),
                mesh = mesh2, spatial = "on")
print(fit2)
## Spatial model fit by ML ['sdmTMB']
## Formula: N ~ Y
## Mesh: mesh2 (isotropic covariance)
## Data: qgl
## Family: poisson(link = 'log')
##  
##             coef.est coef.se
## (Intercept)    -1.90    0.54
## Y               0.24    0.07
## 
## Matérn range: 4.11
## Spatial SD: 0.75
## ML criterion at convergence: 255.395
## 
## See ?tidy.sdmTMB to extract these values as a data frame.

Yの係数は0.24と推定されました。

予測値

予測値を図示します。

p2 <- predict(fit2, newdata = dplyr::select(qgl, X, Y))
qgl$est2 <- exp(p2$est)
ggplot(qgl, aes(x = X, y = Y)) +
  geom_tile(aes(fill = est2)) +
  scale_y_reverse(breaks = c(0, 5, 10)) +
  scale_fill_gradient(name = "株数", low = "grey90", high = "grey10",
                      limits = c(0, 6)) +
  coord_fixed(ratio = 1) +
  theme_bw(base_family = "IPAexGothic")

cutoff = 2とした場合と比較すると、空間的に細かく変動しているようです。