はじめに

令和6年元旦に能登半島で大きな地震がありました．私は能登半島の出身で，そこに実家があるため帰省していました．揺れはとても大きく，街中がぐちゃぐちゃになってしまいました．当然，これだけ大きな被害があると，さまざまな産業に影響が出ると考えられます．農業も例外ではないでしょう．農具が破損したり，農業の担い手が減少したりした地域があると考えられます．　　 この記事では，能登半島の農業集落ごとの水田耕作枚数に注目し，震災前後でそれがどのように増減したかを推定します．
なおこの記事で使用したデータは現在公開しておりません．Sentinel-1とSentinel-2によるリモートセンシングを行い，そのデータを用いてDNNで水田利用状況を推定したものです．もうちょっとしたら公開するかもしれません．

パッケージの読み込み

library(tidyverse) # データハンドリングの効率化
library(sf) # 空間データのハンドリング
library(patchwork) # 複数の図を簡単に並べる
library(agribbit) # 農林業センサスデータを取得
library(ggdist)
library(mapview) # 動かせる地図

データの読み込み

能登半島の水田耕作枚数

predictedという名前が入っている変数は，それぞれの年の水田耕作枚数の予測値．使うのは2021-2024年のデータ．

noto <- read_csv("noto_num_pad_cul.csv") %>% 
  dplyr::select(-...1) # 不要な列を削除
noto %>% glimpse()
## Rows: 806
## Columns: 8
## $ KEY_CODE      <dbl> 1720201001, 1720201002, 1720201004, 1720201005, 17202010…
## $ predicted2017 <dbl> 119, 1, 15, 18, 25, 105, 56, 31, 274, 95, 247, 103, 302,…
## $ predicted2018 <dbl> 146, 1, 19, 23, 39, 112, 56, 29, 268, 91, 264, 114, 302,…
## $ predicted2020 <dbl> 115, 2, 17, 24, 29, 101, 54, 28, 267, 97, 249, 103, 320,…
## $ predicted2021 <dbl> 117, 1, 16, 17, 28, 89, 54, 32, 268, 92, 257, 103, 303, …
## $ predicted2022 <dbl> 120, 0, 18, 17, 36, 105, 58, 28, 283, 88, 254, 109, 299,…
## $ predicted2023 <dbl> 99, 1, 15, 14, 25, 88, 57, 30, 253, 96, 245, 100, 261, 2…
## $ predicted2024 <dbl> 69, 1, 12, 16, 30, 84, 55, 29, 250, 95, 235, 106, 278, 2…

とりあえず描画．{agribbit}で石川県の農林業センサスの.shpを拾います．それを先ほど読み込んだデータと結合して地図にしましょう．agribbit::agri_read_census_shp()の引数に都道府県コード（石川県なら17）を指定すれば取得できます．

ishikawa_shp <- agri.read_census_shp(17) %>% 
  mutate(KEY_CODE = KEY_CODE)

簡単に描画

データを読み込んだら，まずは適当に描画してみるのが鉄の掟，いや，ダイヤモンドの掟です．自分がRを初めて習ったとき（学部2年の冬），その講義の先生がそういう旨のことをおっしゃっていた気がします．
2021年から2024年までの耕作枚数の図のオブジェクトを作成して，それを{patchwork}で並べて表示してみます．

# それぞれの年の図のオブジェクトを作成していく
noto_2021 <- left_join(noto, ishikawa_shp, by = "KEY_CODE") %>% 
  st_as_sf() %>% 
  ggplot() +
  geom_sf(
    aes(fill = predicted2021),
    lwd = .1
  ) +
  scale_fill_fermenter(
    name = "",
    breaks = seq(0, 700, 100),
    direction = 1,
    palette = "YlGnBu"
  ) +
  labs(fill = "", x = "2021") +
  theme_minimal() +
  theme(
    legend.position = "none",
    axis.text = element_blank(),
    axis.ticks = element_blank()
    )

noto_2022 <- left_join(noto, ishikawa_shp, by = "KEY_CODE") %>% 
  st_as_sf() %>% 
  ggplot() +
  geom_sf(
    aes(fill = predicted2022),
    lwd = .1
  ) +
  scale_fill_fermenter(
    name = "",
    breaks = seq(0, 700, 100),
    direction = 1,
    palette = "YlGnBu"
  ) +
  labs(fill = "", x = "2022") +
  theme_minimal() +
  theme(
    legend.position = "none",
    axis.text = element_blank(),
    axis.ticks = element_blank())

noto_2023 <- left_join(noto, ishikawa_shp, by = "KEY_CODE") %>% 
  st_as_sf() %>% 
  ggplot() +
  geom_sf(
    aes(fill = predicted2023),
    lwd = .1
  ) +
  scale_fill_fermenter(
    name = "",
    breaks = seq(0, 700, 100),
    direction = 1,
    palette = "YlGnBu"
  ) +
  labs(fill = "", x = "2023") +
  theme_minimal() +
  theme(
    legend.position = "top",
    legend.key.width = unit(2, "cm"),
    axis.text = element_blank(),
    axis.ticks = element_blank())

noto_2024 <- left_join(noto, ishikawa_shp, by = "KEY_CODE") %>% 
  st_as_sf() %>% 
  ggplot() +
  geom_sf(
    aes(fill = predicted2024),
    lwd = .1
  ) +
  scale_fill_fermenter(
    name = "",
    breaks = seq(0, 700, 100),
    direction = 1,
    palette = "YlGnBu"
  ) +
  labs(fill = "", x = "2024") +
  theme_minimal() +
  theme(
    legend.position = "none",
    axis.text = element_blank(),
    axis.ticks = element_blank())

# 並べる
noto_2021_2024 <- noto_2021 + noto_2022 + noto_2023 + noto_2024 + plot_layout(nrow = 1)
noto_2021_2024

インタラクティブなplotもできる

2024年の水田耕作枚数の地図をインタラクティブに描画してみましょう．任意の場所をクリックすると，その地域の水田耕作枚数を見ることができます．ここで表示される水田耕作枚数は実測値ではなく，あくまでも予測値出ることに注意してください．また，衛星の分解能の問題で小さい水田は正しく観測されていません．そのため．小さい水田は，実際には耕作されていても非耕作水田として分類されている場合があります．この問題を解決するには，より分解能の良い衛星を使うほかありません．分解能の優れた衛星は高額であること，さらに，そのようなデータを持っていてもスピーディに処理できる計算機を有していないことなどの難点がありますがねえ．

noto %>% 
  left_join(., ishikawa_shp) %>% 
  st_as_sf() %>% 
  dplyr::select(predicted2024) %>% 
  mapview(
    zcol = "predicted2024"
  )
## Joining with `by = join_by(KEY_CODE)`

山形県の水田耕作枚数

山形県はDIDにおける統制群として扱う．こちらも2021-2024までのデータしか使わない．

yamagata <- noto_abandoned_df <- read_csv("yamagata_num_pad_cul.csv") %>% 
  mutate(
    KEY_CODE = as.numeric(KEY_CODE)
  ) %>% 
  dplyr::select(-...1)
## New names:
## Rows: 2582 Columns: 7
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," dbl
## (7): ...1, KEY_CODE, predicted2019, predicted2021, predicted2022, predic...
## ℹ 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.
## • `` -> `...1`
yamagata %>% glimpse()
## Rows: 2,582
## Columns: 6
## $ KEY_CODE      <dbl> 620101001, 620101002, 620101003, 620101004, 620101009, 6…
## $ predicted2019 <dbl> 122, 94, 0, 0, 571, 2, 2, 0, 0, 84, 4, 0, 0, 53, 474, 18…
## $ predicted2021 <dbl> 124, 98, 0, 0, 646, 3, 1, 0, 0, 92, 10, 0, 0, 57, 542, 1…
## $ predicted2022 <dbl> 121, 95, 0, 0, 624, 2, 1, 0, 1, 77, 8, 0, 0, 36, 491, 15…
## $ predicted2023 <dbl> 123, 96, 1, 0, 606, 5, 0, 0, 0, 65, 1, 0, 1, 54, 504, 20…
## $ predicted2024 <dbl> 124, 97, 1, 0, 633, 1, 1, 0, 0, 76, 12, 0, 1, 55, 508, 1…

山形も2024年だけ簡単に描画してみるとする．

yamagata_shp <- agri.read_census_shp(6)
yamagata_shp %>% 
  left_join(., yamagata) %>% 
  dplyr::select(predicted2024) %>% 
  mapview(
    zcol = "predicted2024"
  )
## Joining with `by = join_by(KEY_CODE)`

水田耕作枚数の分布

個人的な話ですが，ある変数の値の分布を見る際，カーネル密度推定が好きです．箱ひげ図なども主流ですが，あまり好ましく思っていません．箱ひげずは四分位数や最大値・最小値を拾えますが，分布の形状を教えてくれないからです．分布の形状ってとっても大事な情報だと思うので，それを落としてしまう箱ひげずはちょっと……と感じます．

そんなことはさておき，能登半島の水田耕作枚数の分布を見てみよう．赤い分布が2024年の，背後霊みたいな青い分布が2023年の水田耕作枚数の分布です．こうしてみる震災前後でとあんまり変わっていないような感じがします．しかし，若干2024年の方が，耕作枚数20-30のあたりに大きなボリュームがある気もします．

noto %>% 
  ggplot() +
  geom_density(
    aes(predicted2023),
    fill = "royalblue", color = NA,
    alpha = 1
  ) +
    geom_density(
    aes(predicted2024),
    fill = "tomato", color = NA,
    alpha = .6
  ) +
  labs(x = "the number of cultivatied paddies in Noto Pen.") +
  theme_minimal()

続いて山形県でも見てみましょう．あんまり変化がないですね．山形県でも水田耕作枚数が少ないところのボリュームが，2024年の方で大きくなっているような気がします．

yamagata %>% 
  ggplot() +
  geom_density(
    aes(predicted2023),
    fill = "royalblue", color = NA,
    alpha = 1
  ) +
    geom_density(
    aes(predicted2024),
    fill = "tomato", color = NA,
    alpha = .6
  ) +
  labs(x = "the number of cultivatied paddies in Noto Pen.") +
  theme_minimal()

ちなみにですが，ある巷では，raincloudplotなどというものが支持を得ています．raincloudplotで2021-2024の能登半島の水田耕作枚数の分布を描画してみましょう．

# アドホックにデータフレームをこさえる
hoge1 <- noto$predicted2021
hoge2 <- noto$predicted2022
hoge3 <- noto$predicted2023
hoge4 <- noto$predicted2024
year <- c(rep(2021, length(hoge1)), rep(2022, length(hoge1)), rep(2023, length(hoge1)), rep(2024, length(hoge1)))
hoge <- data_frame(
  val = c(hoge1, hoge2, hoge3, hoge4),
  year = year
)

# 描画
hoge %>% 
  ggplot() +
  stat_halfeye(
    aes(x = val, y = as_factor(year),
        fill = as_factor(year))
  ) +
  geom_boxplot(
    aes(x = val, y = as_factor(year), fill = as_factor(year)),
    position = position_nudge(y = -.1),
    width = .1,
    outlier.colour = NA,
    alpha = .5
  ) +
  geom_point(
    aes(x = val, y = as_factor(year)),
    color = "black",
    position = position_jitter(width = 0, height = 0.025, seed = 1),
    alpha = .05
  ) +
  xlim(0, 400) +
  scale_fill_manual(
    values = c("tomato", "royalblue", "yellow", "pink2")
  ) +
  labs(x = "the numner of paddies cultivarted", y = "") +
  theme_minimal() +
  theme(legend.position = "none")

raincloudplotとは，こんな感じのやつです↑ 『ぼっち・ざ・ろっく』を意識したカラーリングをしてみました．このような描画方法では，分布の形状と四分位数の情報を同時に表現できますね．いいね！
描画はこんくらいで終わっておきましょう．

DIDやで

DID概説

さて，令和6年能登半島地震によって，能登半島の水田耕作枚数がどの程度減少したのかを推定してみましょう．DIDをやります．
DIDが何かということについては他の文献に譲ります．が，簡単に説明すると，以下の式の係数\(\beta_3\)を知ろうとするものです．\(time_i\)と\(treat_i\)はそれぞれダミー変数です．\(time_i\)は地域\(i\)がの観測値が2024年のものなら1，2023年のものなら0をとります．\(treat_i\)は介入の有無を表しており，iが能登半島の農業集落であれば1，山形県の集落であれば0を取ります．
\[ y_i = \beta_0 + \beta_1time_i + \beta_2tread_i + \beta_3time_i treat_i \] それぞれの係数の解釈に関しては，以下の図をご覧ください．ざっくりいうと\(\beta_3\)は，震災がなかった世界線における2024年の能登半島の水田耕作枚数（これは観測できないが，山形県のデータを元にして再現）と，震災があった世界線における能登半島の水田耕作枚数（これは観測できる）の差です．

DID実装

統制群の選抜

さて，実際にDIDを行います．能登半島の各集落(806地域)の震災の影響を推定したいので，806回のDIDを行います．
DIDでは平行トレンド仮定が満たされていることが望ましいです．そこで，震災が発生するまで（2021-2023）の水田耕作枚数の傾向が似ている山形県の集落を選抜して，統制群として使用します．ある能登半島の集落\(i\)に対して，以下の式を小さくする山形県の上位10地域を，統制群とします．ただし\(p_{it}\)は集落\(i\)の\(t\)年の水田耕作枚数を表します． \[ diff_{ij} = \sum_{t=2021}^{2023}(p_{it} - p_{jt})^2 \]

このような操作をすることで，比較対象として望ましい地域を用いることができると思われます．すなわち，平行トレンド仮定が満たされるのではないかと思います．選抜の例を出します．青線は能登半島のある地域の水田耕作枚数の推移です．細い緑線は選抜された山形県の集落10地域の推移です．太い緑線は，山形県10地域の平均値です．消え入るような黒い点線は，選抜されなかった山形県の集落たちの推移です．

806回のDID

effect_vector <- c()
KEY_CODE <- c()
for (i in 1:nrow(noto)){
  noto_extracted <- noto[i,] # 検証したい能登半島の地域を抜き出す
  treat_2023_value <- noto_extracted$predicted2023 # 基準となる地域の耕作枚数2023
  treat_2022_value <- noto_extracted$predicted2022 # 基準となる地域の耕作枚数2022
  treat_2021_value <- noto_extracted$predicted2021 # 基準となる地域の耕作枚数2021
  KEY_CODE <- c(KEY_CODE, noto_extracted$KEY_CODE)
  # ここでyamagataについてのjループ
  # 基準地域との耕作枚数のmseが小さい10地域を選抜
  yamagata_sum_vec <- c()
  for (j in 1:nrow(yamagata)){
    mse21 <- (treat_2021_value - yamagata[j,]$predicted2021)^2
    mse22 <- (treat_2022_value - yamagata[j,]$predicted2022)^2
    mse23 <- (treat_2023_value - yamagata[j,]$predicted2023)^2
    mse_sum <- mse21 + mse22 + mse23
    yamagata_sum_vec <- c(yamagata_sum_vec, mse_sum)
  }
  yamagata_key_code <- yamagata$KEY_CODE %>% as.vector()
  thereshold <- sort(yamagata_sum_vec)[10] # mseが小さい上位10地域を選抜
  df <- data.frame(yamagata_key_code, yamagata_sum_vec) %>% 
    filter(yamagata_sum_vec <= thereshold) %>% 
    rename(KEY_CODE = yamagata_key_code)
  # ここで作成したdfを元にして，parallelを満たすyamagataのデータフレームをこさえる
  yamagata_extracted <- left_join(df, yamagata, by = "KEY_CODE")
  # ここからDIDの本編
  if (nrow(yamagata_extracted) != 0){
    # ここから選択したデータの整理を行う
    # 従属変数
    noto_cult_2023 <- c(noto_extracted$predicted2023)
    noto_cult_2024 <- c(noto_extracted$predicted2024)
    yamagata_cult_2023 <- c(yamagata_extracted$predicted2023)
    yamagata_cult_2024 <- c(yamagata_extracted$predicted2024)
    y <- c(
      noto_cult_2023,
      noto_cult_2024,
      yamagata_cult_2023,
      yamagata_cult_2024
    )
    # 説明変数
    # 介入ダミー
    treat  <- c(
      rep(1, length(noto_cult_2023)),
      rep(1, length(noto_cult_2024)),
      rep(0, length(yamagata_cult_2023)),
      rep(0, length(yamagata_cult_2023))
    )
    # 時間ダミー
    time <- c(
      rep(0, length(noto_cult_2023)),
      rep(1, length(noto_cult_2024)),
      rep(0, length(yamagata_cult_2023)),
      rep(1, length(yamagata_cult_2023))
    )
    # データフレーム化
    df <- data.frame(
      y,
      treat,
      time
    )
    # ここから回帰
    formula <- y ~ treat*time # 回帰式の定義
    did_model <- lm(formula, df)
    effect_vector <-c(effect_vector, as.numeric(did_model$coefficients[4]))
  }
  else{
    effect_vector <-c(effect_vector, NA)
  }
  if (i %% 100 == 0){
    print("processing")
  }
}
## [1] "processing"
## [1] "processing"
## [1] "processing"
## [1] "processing"
## [1] "processing"
## [1] "processing"
## [1] "processing"
## [1] "processing"
# 新型DIDでの結果チェック
did_result_df <- data.frame(
  KEY_CODE,
  effect_vector
)

ちょっと時間がかかるかもしれませんが，これで能登半島806地域のDIDが完了しました．

DIDの結果

最後にDIDの結果を見てみましょう！ 地図にしてみます．

left_join(noto, did_result_df, by = "KEY_CODE") %>% 
  left_join(., ishikawa_shp) %>% 
  st_as_sf() %>% 
  st_simplify(dTolerance = 5) %>% 
  ggplot() +
  geom_sf(
    aes(fill = effect_vector),
    lwd = .1) +
  scale_fill_gradient2(
    low = "tomato",
    mid = "white", 
    high = "seagreen3"
  ) +
  labs(fill = "DID effect") +
  theme_minimal()
## Joining with `by = join_by(KEY_CODE)`

赤色の地域は，震災で耕作枚数が減少したと考えられる地域です．緑色の地域は耕作枚数が増えたと考えられる地域です． この結果を見た感じ，北部と南部で傾向に違いがあることがわかります．北部は減少傾向にあるのに対して，南部では増加傾向になります．実際，震災の影響が大きかったのは北部です．そう考えると，直感に妥当する結果ではあります．しかし，増えた地域があるというのは自分的には懐疑的です．

さいごに

能登半島が，また，そこに住むみんなが，幸せに暮らせる場所になればいいと思います．自分の研究はそれにはあまり貢献できませんが，遠くから応援したいと思います．