sfdep で手軽に空間統計

# sfdep で手軽に空間統計

R の `sfdep` パッケージを使うと割と手軽に空間統計できることに気付いたのでメモ。
https://cran.r-project.org/web/packages/sfdep/ 

パッケージの準備はこんな感じ。

```r
if(!require("tidyverse")) {
    install.packages(c("tidyverse"))
    require(tidyverse)
}

if(!require("sfdep")) {
    install.packages(c("sfdep"))
    require(sfdep)
}

if(!require("NipponMap")) {
    install.packages(c("NipponMap"))
}

library(tidyverse)
```

NipponMap から日本のシェープファイル（人口の情報付き）を抜き出す。

```r
shapefile <- 
    "shapes/jpn.shp" |>
    system.file(package = "NipponMap") |>
    sf::read_sf() |> 
    sf::st_set_crs(4326)
```

変数を設定しておく。
※sfdep の公式の解説 https://sfdep.josiahparry.com/ では sfdep::st_knn ではなくsfdep::st_contiguity を扱っている。題材がフランスだからそれで成立しているけど、日本の都道府県でそれをやると北海道と沖縄が他の都道府県と「隣接」していない関係でエラーが出るので、今回は sfdep::st_knn を使用しておく。

```r
nb_func = {\(x) x |> sf::st_centroid() |> sfdep::st_knn(k = 4)}  # 隣接行列
nsim = 1e5  # Local Moran の試行回数（再現性のため増やしておく）
p_threshold = 0.05  # Local Moran の有意水準

factor_levels = c(
    "High-High", "High-Low", "Low-High", "Low-Low",
    "Not Significant"
)

factor_name = "Quadrant"

label_style = scales::label_number(scale_cut = scales::cut_short_scale())

# 散布図の設定

scatter_color_palette = c(
    `High-High` = "red", `High-Low` = "salmon", 
    `Low-High` = "lightblue", `Low-Low` = "blue",
    `Not Significant` = "black"
)

scatter_shape_palette = c(
    `High-High` = 20, `High-Low` = 20, `Low-High` = 20, `Low-Low` = 20,
    `Not Significant` = 1
)

scatter_size_palette = c(
    `High-High` = 4, `High-Low` = 4, `Low-High` = 4, `Low-Low` = 4,
    `Not Significant` = 2
)

x_label = "Population"
y_label = "Spatial Lag"

# 地図の設定

crs_value = 6677

lisa_color_palette = c(
    `High-High` = "red", `High-Low` = "salmon", 
    `Low-High` = "lightblue", `Low-Low` = "blue",
    `Not Significant` = "white"
)
```

データの下ごしらえはこんな感じ。

```r
jpn_map <-
    shapefile |> 
    dplyr::mutate(
        nb = geometry |> nb_func(),
        wt = nb |> sfdep::st_weights(),
        spatial_lag = population |> sfdep::st_lag(nb, wt),
        moran = population |> sfdep::local_moran(nb, wt, nsim = nsim),
    ) |> 
    tidyr::unnest(moran) |> 
    dplyr::mutate(
        pysal = dplyr::if_else(
            p_folded_sim < p_threshold,
            as.character(pysal),
            "Not Significant"
        ) |> 
            factor(levels = factor_levels)
    )
```

Local Moran *I* の散布図ならこんな感じ。

```r
p_scatter <- 
    jpn_map |> 
    ggplot2::ggplot(aes(x = population, y = spatial_lag)) +
    geom_point(aes(color = pysal, shape = pysal, size = pysal)) +
    geom_vline(aes(xintercept = mean(population)), lty = 2) + 
    geom_hline(aes(yintercept = mean(spatial_lag)), lty = 2) + 
    scale_x_continuous(labels = label_style) +
    scale_y_continuous(labels = label_style) +
    scale_color_manual(
        name = factor_name,
        values = scatter_color_palette,
    ) +
    scale_shape_manual(
        name = factor_name,
        values = scatter_shape_palette,
    ) +
    scale_size_manual(
        name = factor_name,
        values = scatter_size_palette,
    ) +
    labs(x = x_label, y = y_label) +
    coord_fixed(ylim = c(min(jpn_map$population), max(jpn_map$population))) +
    cowplot::theme_cowplot() +
    cowplot::background_grid()
```
![scatter_plot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/346391/5a87dbc9-4892-0136-f021-43b33a447885.png)


人口分布ならこんな感じ。

```r
p_population <- 
    jpn_map |> 
    ggplot2::ggplot() +
    geom_sf(aes(fill = population)) +
    coord_sf(crs = crs_value) +
    scale_fill_viridis_c(name = "Population", labels = label_style) + 
    cowplot::theme_cowplot() +
    cowplot::background_grid()
```
![population_map.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/346391/24d703f2-88b4-a4a1-ea98-f68eb3985de3.png)


LISA cluster ならこんな感じ。

```r
p_lisa <- 
    jpn_map |> 
    ggplot2::ggplot() +
    geom_sf(aes(fill = pysal)) +
    coord_sf(crs = crs_value) +
    scale_fill_manual(
        name = factor_name,
        values = lisa_color_palette,
    ) +
    cowplot::theme_cowplot() +
    cowplot::background_grid()
```
![lisa_plot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/346391/bf2b1324-74af-a06a-e51c-50f3b8171863.png)