15 生態学

必須パッケージ

この章では、Chapter 2 ~ Chapter 5 で説明した地理データ解析と処理について十分に理解していることを前提に説明する。 また、、Chapters 10 で解説した R の専用 GIS インタフェースや、Chapter 12 の空間クロスバリデーション(CV)も活用することができる。

この章では、以下のパッケージを使用する。

library(sf)
library(terra)
library(dplyr)
library(data.table)        # fast data.frame manipulation (used by mlr3)
library(mlr3)              # machine learning (see Chapter 12)
library(mlr3spatiotempcv)  # spatio-temporal resampling 
library(mlr3tuning)        # hyperparameter tuning package
library(mlr3learners)      # interface to most important machine learning packages
library(paradox)           # defining hyperparameter spaces
library(ranger)            # random forest package
library(qgisprocess)       # bridge to QGIS (Chapter 10)
library(tree)              # decision tree package
library(vegan)             # community ecology package

15.1 イントロダクション

本章では、霧のオアシスの植生勾配をモデル化し、水の利用可能性に明らかに支配されている特徴的な植生帯を明らかにする。 そのために、これまでの章で紹介した概念をまとめ、さらに拡張していく(Chapter 25 、Chapters 10、Chapter 12)。

霧のオアシスは、これまでに出会った中で最も魅力的な植生形態の一つである。 このような地層は、地元では lomas と呼ばれ、ペルーやチリの海岸砂漠に沿った山々に発達している。96 砂漠の極限状態と遠隔地であることが、霧のオアシスの固有種を含むユニークな生態系の生息地となっている。 年間平均 30〜50 mm 程度の降水量と乾燥した条件にもかかわらず、霧の発生により、南半球の冬の間に植物が利用できる水量が増加する。 その結果、ペルーの海岸線に沿った南向きの山の斜面が緑色になる(Figure 15.1)。 冬のフンボルト海流がもたらす逆転現象の下に発生する霧が、この生息地の名前の由来となっている。 数年に一度、エルニーニョ現象によって、この太陽の降り注ぐ環境に集中豪雨がもたらされる (Dillon, Nakazawa, and Leiva 2003)。 その結果、砂漠に花が咲き、木の苗が次の乾燥した環境を生き抜くために十分な長さの根を張ることが可能になる。

残念ながら、霧のオアシスは、主に人間の活動(農業と気候変動)により、大きく危機に瀕している。 このユニークな植生生態系の最後の残りを効果的に保護するためには、原生植物相の構成と空間分布に関するエビデンスが必要である (Muenchow, Bräuning, et al. 2013; Muenchow, Hauenstein, et al. 2013)。 また、 lomas のある山には観光地としての経済的な価値もあり、レクリエーションを通じて地域住民の幸福にも貢献する。 例えば、ペルー人の多くは海岸沿いの砂漠に住んでいて、「緑」に一番近いのはロマス山であることが多いんである。

この章では、前の章で学んだ技術を生態学的に実際に応用していこう。 事例として、ペルー中央北岸のカスマ近郊にある lomas のある山、モンゴン山の南斜面における維管束植物の構成と空間分布を分析する(Figure 15.1)。

Mongón 山調査地、Muenchow, Schratz, and Brenning (2017) より。

FIGURE 15.1: Mongón 山調査地、Muenchow, Schratz, and Brenning (2017) より。

モンゴル山への野外調査において、2011年の冬、ランダムにサンプルした 4x4 m2 の 100 区画に生息するすべての維管束植物を記録した (Muenchow, Bräuning, et al. 2013)。 サンプル採取は、その年の強いラニーニャ現象に重なった(ENSO monitoring of NOASS Climate Prediction Centerを参照)。 そのため、沿岸部の砂漠は通常よりもさらに高い乾燥度を示すようになった。 一方、ペルーの lomas のある山脈の南斜面では、霧の活動も活発になっている。

Ordinations は、(ノイズの多い)データセットから主要な勾配を抽出するための次元削減技術であり、私たちの場合は南山斜面に沿って発達した植物学的勾配である(次のセクションを参照)。 本章では、最初の序列軸である植物相の勾配を、標高、傾斜、集水域および NDVI といった環境予測因子の関数としてモデル化することにする。 このために、非常に有名な機械学習アルゴリズムであるランダムフォレストモデルを利用する (Breiman 2001)。 このモデルを使えば、調査地のどこでも植物組成の空間的な予測が可能になる。 最適な予測を保証するために、空間クロスバリデーション (Section 12.5.2 参照)の助けを借りて、ハイパーパラメータを事前に調整することが推奨される。

15.2 データとデータ準備

以降の解析に必要なデータは、spDataLargeパッケージから入手可能である。

data("study_area", "random_points", "comm", package = "spDataLarge")
dem = rast(system.file("raster/dem.tif", package = "spDataLarge"))
ndvi = rast(system.file("raster/ndvi.tif", package = "spDataLarge"))

study_area は調査地域の外形を表すポリゴン、 はランダムに選ばれた 100 地点を含む オブジェクトである。 random_points sf comm はワイドデータ形式 の群集行列で、行はフィールドで訪問した場所、列は観察された種を表している。97

# 35から40のサイトと、それに対応する群集マトリックスの
# 最初の5種の発生状況
comm[35:40, 1:5]
#>    Alon_meri Alst_line Alte_hali Alte_porr Anth_eccr
#> 35         0         0         0       0.0     1.000
#> 36         0         0         1       0.0     0.500
#> 37         0         0         0       0.0     0.125
#> 38         0         0         0       0.0     3.000
#> 39         0         0         0       0.0     2.000
#> 40         0         0         0       0.2     0.125

数値はサイトごとの種の被度を表し、サイト面積に対する種の被度の割合(%;ひとつのサイトでは、個々の植物間の被度の重複により100%を超える場合があることに注意)で記録された。 comm の rownames は random_pointsid 列に対応する。 dem は調査地域の数値標高モデル (DEM)は Landsat シーンの赤色および近赤外チャネルから算出した正規化植生指標(NDVIndvi)である(Section 4.3.3 および ?spDataLarge::ndvi.tif を参照)。 Figure 15.2 に示すように、demrandom_pointsstudy_area を重ねて表示すことで、データをより身近なものにすることができる。

スタディマスク(ポリゴン)、サンプリング地点(黒点)、背景の DEM。

FIGURE 15.2: スタディマスク(ポリゴン)、サンプリング地点(黒点)、背景の DEM。

次のステップは、モデリングと予測地図作成に必要な変数( Section 15.4.2 参照)だけでなく、Non-metric multidimensional scaling (NMDS) 軸を調査地域の主要勾配、高度と湿度にそれぞれ整合させるための変数も計算することである(Section 15.3 参照)。

具体的には、R-GIS ブリッジを用いて、デジタル標高モデルから集水勾配と集水面積を計算する(Chapter 10 参照)。 曲率も重要な予測因子である可能性があり、「演習」セクションで、それらがモデリング結果にどのような影響を与えるかを確認することができる。

集水域と集水勾配を計算するには、saga:sagawetnessindex 関数を利用することができる。98 qgis_show_help() 特定のジオアルゴリズムのすべての関数パラメータとデフォルト値を返す。 ここでは、その一部のみを紹介する。

qgisprocess::qgis_show_help("saga:sagawetnessindex")
#> Saga wetness index (saga:sagawetnessindex)
#> ...
#> ----------------
#> Arguments
#> ----------------
#> 
#> DEM: Elevation
#>  Argument type:  raster
#>  Acceptable values:
#>      - Path to a raster layer
#> ...
#> SLOPE_TYPE: Type of Slope
#>  Argument type:  enum
#>  Available values:
#>      - 0: [0] local slope
#>      - 1: [1] catchment slope
#> ...
#> AREA: Catchment area
#>  Argument type:  rasterDestination
#>  Acceptable values:
#>      - Path for new raster layer
#>... 
#> ----------------
#> Outputs
#> ----------------
#> 
#> AREA: <outputRaster>
#>  Catchment area
#> SLOPE: <outputRaster>
#>  Catchment slope
#> ...

次に、R の名前付き引数を使って、必要なパラメータを指定する(Section 10.2 を参照)。 R のグローバル環境に住んでいる SpatRaster を使って、入力ラスタ DEM を指定できることを思い出そう(Section 10.2 参照)。 SLOPE_TYPE に 1 を指定することで、アルゴリズムが集水勾配を返すようになる。 出来上がったラスタ は、SAGA のラスタフォーマットである .sdat という拡張子で一時ファイルに保存される。

# 環境予測因子:集水勾配と集水面積
ep = qgisprocess::qgis_run_algorithm(
  alg = "saga:sagawetnessindex",
  DEM = dem,
  SLOPE_TYPE = 1, 
  SLOPE = tempfile(fileext = ".sdat"),
  AREA = tempfile(fileext = ".sdat"),
  .quiet = TRUE)

これは、計算された出力ラスタへのパスを含む ep という名前のリストを返す。 集水域と集水勾配を多層構造 SpatRaster オブジェクトに読み込んでみよう(Section 2.3.4 参照)。 さらに、これに2つのラスタオブジェクト、すなわち demndvi を追加する。

# 集水域と集水勾配を読み取る
ep = ep[c("AREA", "SLOPE")] |>
  unlist() |>
  terra::rast()
names(ep) = c("carea", "cslope") # 名前をつける
terra::origin(ep) = terra::origin(dem) # ラスタが同じ原点を持つことを確認
ep = c(dem, ndvi, ep) # 多層 SpatRaster オブジェクトに dem と ndvi を追加する。

さらに、集水域 の値は右側に大きく偏っている (hist(ep$carea))。 対数10変換をすると、より正規分布に近くなる。

ep$carea = log10(ep$carea)

読者の便宜を図るため、spDataLargeep を追加した。

ep = terra::rast(system.file("raster/ep.tif", package = "spDataLarge"))

最後に、現地観測に地形属性を抽出することができる(Section 6.3 も参照)。

# terra::extract は、不要な ID 列を自動的に追加
ep_rp = terra::extract(ep, random_points) |>
  dplyr::select(-ID)
random_points = cbind(random_points, ep_rp)

15.3 次元性を低減

順序付け(Ordination)は、植生学において、0 で埋め尽くされた大規模な種間プロット行列から主要情報(生態的勾配に相当することが多い)を抽出するための一般的なツールである。 しかし、リモートセンシング、土壌学、ジオマーケティング などの分野でも利用されている。 順序付けテクニックに馴染みがない場合、または復習が必要な場合は、生態学で人気の順序付けテクニックを簡単に紹介した Michael W. Palmer の web page と R でこれらのテクニックを適用する方法について深く調べた Borcard, Gillet, and Legendre (2011) に目を通してみよう。 vegan のパッケージのドキュメントも非常に有用なリソースである (vignette(package = "vegan"))。

主成分分析 (principal component analysis, PCA) は、おそらく最も有名な順序付け の手法である。 変数間の線形関係が期待でき、2 つのプロット(オブザベーション)における変数の共同不在が類似性とみなせる場合、次元を削減するためのすばらしいツールである。 これは植生データではほとんどない。

ひとつは、植物の存在は通常、勾配(湿度、温度、塩分など)に沿って、最も好ましい条件でピークを迎え、好ましくない条件に向かって減少していくという単峰性の関係にある。

第二に、ある種が 2 つの区画で同時に存在しないことは、類似性を示す指標とはなりにくい。 ある植物種が、サンプルの中で最も乾燥した場所(例:極度の砂漠)と最も湿った場所(例:木のサバンナ)から姿を消したとする。 というのも、この2つの全く異なる環境設定に共通するのは、(稀少なユビキタス種を除いて)種が存在しないという点だけである可能性が非常に高いからである。

Non-metric multidimensional scaling (NMDS) は、生態学でよく使われる次元削減手法の一つである (von Wehrden et al. 2009)。 NMDS は、元の行列のオブジェクト間の距離と、順序付けられたオブジェクト間の距離の間のランクベースの差異を低減する。 その差をストレスとして表現している。 ストレス値が低いほど、順序付け、すなわち元の行列の低次元表現が良好であることを示す。 ストレス値が 10 より小さいと適合性が高く、15程度でも良好、20より大きいと適合性が低いことを表している (McCune, Grace, and Urban 2002)。 R では、vegan パッケージの metaMDS() で NMDS を実行することができる。 入力として、サイトを行、種を列とする群集行列を期待する。 しばしば、有-無データを用いた順序付け は、(説明される分散の点で)より良い結果をもたらするが、その代償として、もちろん、入力行列の情報量は少なくなる(演習も参照)。 decostand() は、数値の観測結果を、1が種の発生、0が種の不在を示す有無に変換する。 NMDS のようなオーダリング技術では、部位ごとに少なくとも1回の観測が必要である。 したがって、種が発見されなかったサイトはすべて除外する必要がある。

# presence-absence 行列
pa = vegan::decostand(comm, "pa") # 100 行(箇所), 69 列(種)
# 1種以上発見されたサイトのみを残す
pa = pa[rowSums(pa) != 0, ]  # 84 行, 69 列(種)

得られた行列は、NMDS の入力として機能する。 k は出力軸数を指定し、ここでは 4 とする。99 NMDS は、各ステップで順序付けされた空間をより入力行列に近づけようとする反復的な手順である。 アルゴリズムの収束を確認するために、ステップ数を 500 に設定した(try パラメータ)。

set.seed(25072018)
nmds = vegan::metaMDS(comm = pa, k = 4, try = 500)
nmds$stress
#> ...
#> Run 498 stress 0.08834745 
#> ... Procrustes: rmse 0.004100446  max resid 0.03041186 
#> Run 499 stress 0.08874805 
#> ... Procrustes: rmse 0.01822361  max resid 0.08054538 
#> Run 500 stress 0.08863627 
#> ... Procrustes: rmse 0.01421176  max resid 0.04985418 
#> *** Solution reached
#> 0.08831395

ストレス値 9 は非常に良い結果を表し、縮小された順序付け空間が入力行列の分散の大部分を表していることを意味する。 全体として、NMDS は、順序付け空間において、(種の構成という点で)より類似したオブジェクトがより近くに配置される。 しかし、他の多くの順序付け の手法とは対照的に、軸は任意であり、必ずしも重要度 (Borcard, Gillet, and Legendre 2011) によって順序付けされるわけではない。 しかし、調査地では湿度が主な勾配を表していることが既に分かっている (Muenchow, Bräuning, et al. 2013; Muenchow, Schratz, and Brenning 2017)。 湿度は標高と高い相関があるため、標高に応じて NMDS 軸 を回転させる(NMDS 軸の回転の詳細については ?MDSrotate も参照)。 結果をプロットしてみると、意図したとおり、第1軸は明らかに高度(Figure 15.3)と関連していることがわかる。

elev = dplyr::filter(random_points, id %in% rownames(pa)) |> 
  dplyr::pull(dem)
# 高度(湿度の代理)に応じて NMDS を回転
rotnmds = vegan::MDSrotate(nmds, elev)
# 最初の2軸の抽出
sc = vegan::scores(rotnmds, choices = 1:2, display = "sites")
# 第1軸を高度にプロット
plot(y = sc[, 1], x = elev, xlab = "elevation in m", 
     ylab = "First NMDS axis", cex.lab = 0.8, cex.axis = 0.8)
NMDS の第1軸を高度に対してプロット。

FIGURE 15.3: NMDS の第1軸を高度に対してプロット。

最初の NMDS 軸のスコアは、Mont.Mongón の斜面に沿って現れる異なる植生形態、すなわち植物学的勾配を表している。 それらを空間的に可視化するために、NMDS のスコアを先に作成した予測因子(Section 15.2)でモデル化し、得られたモデルを予測マッピングに利用する(次節参照)ことができる。

15.4 植物相の勾配をモデル化

植物相の勾配を空間的に予測するために、ランダムフォレスト モデル (Tomislav Hengl et al. 2018) を使用する。 ランダムフォレスト モデルは、環境・生態系のモデリングに頻繁に適用され、予測性能の面で最良の結果をもたらすことが多い (Schratz et al. 2019)。 ここで、決定木とバギングについて簡単に紹介する。これらはランダムフォレストの基礎を形成するものであるため、 。 ランダムフォレスト と関連する技術についてのより詳細な説明は James et al. (2013) を参照。

決定木を例として紹介すると、まず、回転した NMDS スコアと現場観測値 (random_points) を結合して、応答予測行列を構築する。 また、得られたデータフレームは、後で mlr3 のモデリングに使用する予定である。

# response-predictor 行列を作成
# id- と response 変数
rp = data.frame(id = as.numeric(rownames(sc)), sc = sc[, 1])
# 予測因子(dem、ndvi、terrain 属性)を結合 
rp = inner_join(random_points, rp, by = "id")

決定木は、予測変数空間をいくつかの領域に分割する。 これを説明するために、最初の NMDS の軸のスコアを応答(sc)、高度(dem)を唯一の予測因子として、このデータに決定木を適用してみる。

tree_mo = tree::tree(sc ~ dem, data = rp)
plot(tree_mo)
text(tree_mo, pretty = 0)
3つの内部ノードと4つの終端ノードを持つ決定木の単純な例。

FIGURE 15.4: 3つの内部ノードと4つの終端ノードを持つ決定木の単純な例。

結果として得られる木は、3つの内部ノードと4つの終端ノード( Figure 15.4 )で構成される。 木の一番上にある最初の内部ノードは、328.5 m 以下のすべてのオブザベーションを左に、それ以外のすべてのオブザベーションを右の枝に割り当てる。

左の枝に入るオブザベーションは、平均 NMDS スコアが -1.198 である。 -1.198. 全体として、標高が高いほどNMDS のスコアが高くなる、というように解釈できる。 つまり、単純な決定木によって、すでに4つの異なる植物群像が明らかにされているのである。 詳しく学びたい方は、Section 15.4.2 を参照。 決定木は に過剰に適合する傾向がある。つまり、ノイズを含む入力データを忠実に反映しすぎるため、予測性能が低下するのである。 [ Section 12.4 ; James et al. (2013)] . ブートストラップ集計(バギング)は、この問題を克服するためのアンサンブル手法である。 アンサンブル技術は、複数のモデルの予測値を単純に結合するものである。 このように、バギングでは同じ入力データから繰り返しサンプルを取り、その予測値を平均化する。 これにより、分散とオーバーフィッティング、決定木と比較してはるかに優れた予測精度を実現する。 最後に、ランダムフォレストは、相関の高い木の予測を平均化すると、相関の低い木の予測を平均化するよりも分散が大きく、信頼性が低くなるので、バギングを拡張して改良することが望ましい (James et al. 2013)。 これを実現するために、ランダムフォレストはバギングを使用するが、従来のバギングでは各木が利用可能なすべての予測子を使用できるのとは対照的に、ランダムフォレストは利用可能なすべての予測子のランダムサンプルだけを使用する。

15.4.1 mlr3 のビルドブロック

このセクションのコードは、Section 12.5.2 で紹介したステップをほぼ踏襲している。 違いは以下の通りである。

  1. 応答変数は数値なので、回帰 タスクは、Section 12.5.2 の分類 タスクに取って代わる。
  2. カテゴリー応答変数にしか使えないAUROC の代わりに、性能指標として平均二乗誤差(RMSE )を使うことにする。
  3. サポートベクタマシン の代わりに、ランダムフォレスト モデルを使用している。これは当然ながら、異なるハイパーパラメータを伴う。
  4. バイアスを低減した性能指標の評価は、読者の皆様への課題として残している(演習問題参照)。 その代わりに、(空間)予測のためのハイパーパラメータを調整する方法を示す。

100 回繰り返しの 5 回空間クロスバリデーション と 50 回のランダムサーチを使用した場合、バイアス低減された性能推定値を得るために 125,500 個のモデルが必要だったことを思い出してみよう (Section 12.5.2 を参照)。 ハイパーパラメータ のチューニングレベルでは、最適なハイパーパラメータの組み合わせを見つけ、それを特定の空間分割のテストデータを予測するための外部パフォーマンスレベルで使用した(Figure 12.6 も参照)。 これを 5 つの空間分割に対して行い、100 回繰り返した結果、合計 500 の最適なハイパーパラメータの組み合わせが得られた。 空間予測に使うべきはどちらか? 答えは簡単で、全くない。 このチューニングは、バイアスを低減した性能推定値を得るために行われたものであり、最良の空間予測を行うために行われたものではないことに留意してみよう。 後者については、完全なデータセットから最適なハイパーパラメータの組み合わせを推定する。 これは、内部のハイパーパラメータ のチューニング・レベルがもはや必要ないことを意味する。これは、真のアウトカムが利用できない新しいデータ(未訪問のフィールド観測)に私たちのモデルを適用するので、完全に理にかなっており、したがってテストはいかなる場合でも不可能なのである。 そこで、5回繰り返しの空間 CV によって、完全なデータセットで良好な空間予測を行うためにハイパーパラメータを調整することにした。

既に入力変数 (rp) を構築しているので、mlr3 の構成要素(タスク、学習器、リサンプリング)を指定するための準備は全て整っている。 空間タスクの指定には、再び mlr3spatiotempcv パッケージを使用する (Schratz et al. 2021 & Section 12.5)。そして、私たちの応答 (sc) は数値なので、回帰タスクを使用する。

# task を作成
task = mlr3spatiotempcv::as_task_regr_st(dplyr::select(rp, -id, -spri),
  id = "mongon", target = "sc")

バックエンドとして sf オブジェクトを使用すると、後の空間分割に必要なジオメトリ情報が自動的に提供される。 さらに、idspri の列は、これらの変数をモデリングにおける予測因子として使用すべきではないため、削除した。 次に、ranger パッケージのランダムフォレスト 学習器を構築する (Wright and Ziegler 2017)

lrn_rf = lrn("regr.ranger", predict_type = "response")

例えばサポートベクタマシン(Section 12.5.2 参照)とは対照的に、ランダムフォレストはハイパーパラメータのデフォルト値で使用した場合、既に良い性能を示すことが多い(これが人気の理由の一つだろう)。 それでも、チューニングによってモデル結果が適度に改善されることが多いので、努力する価値はある (Probst, Wright, and Boulesteix 2018)。 ランダムフォレストでは、ハイパーパラメータ mtrymin.node.sizesample.fraction がランダム性の度合いを決定するので、これらを調整する必要がある (Probst, Wright, and Boulesteix 2018)mtry は、各ツリーでいくつの予測変数を使用すべきかを示す。 すべての予測変数が使用される場合、これは事実上バギングに相当する(Section 15.4 の冒頭を参照)。 sample.fraction パラメータは、各ツリーで使用されるオブザベーションの割合を指定する。 分画が小さいと多様性が増すので、相関のある樹木が少なくなり、望ましいことが多い(上記参照)。 min.node.size パラメータは、端末ノードが少なくとも持つべき観測値の数を示す( Figure 15.4 も参照)。 当然ながら、木や演算時間が大きくなればなるほど、min.node.size は小さくなる。

ハイパーパラメータの組み合わせはランダムに選択されるが、特定のチューニング限界(paradox::ps() で作成)の範囲内に収まる必要がある。 mtry は 1 から予測変数の数 4) までの範囲でなければならない。sample.fraction は 0.2 から 0.9 の間、min.node.size は 1 から 10 の間でなければならない (Probst, Wright, and Boulesteix 2018)

# 探索空間を指定
search_space = paradox::ps(
  mtry = paradox::p_int(lower = 1, upper = ncol(task$data()) - 1),
  sample.fraction = paradox::p_dbl(lower = 0.2, upper = 0.9),
  min.node.size = paradox::p_int(lower = 1, upper = 10)
)

探索空間を定義したことで、AutoTuner() 関数でチューニングを指定する準備が整いた。 地理的なデータを扱うので、今回も空間クロスバリデーションを用いてハイパーパラメータを調整する(Section 12.4 と Section 12.5 を参照)。 具体的には、1回だけ繰り返す5分割の空間分割(rsmp())を使用することにする。 これらの空間分割のそれぞれにおいて、あらかじめ定義された限界値(seach_space)の範囲内でランダムに選択されたハイパーパラメータ構成( tnr() )を用いながら50個のモデル(trm())を実行し、最適なハイパーパラメータの組合せを見出す (Section 12.5.2https://mlr3book.mlr-org.com/optimization.html#autotuner, Becker et al. 2022 を参照)。 性能指標は二乗平均平方根誤差(RMSE)である。

autotuner_rf = mlr3tuning::AutoTuner$new(
  learner = lrn_rf,
  resampling = mlr3::rsmp("spcv_coords", folds = 5), # spatial partitioning
  measure = mlr3::msr("regr.rmse"), # performance measure
  terminator = mlr3tuning::trm("evals", n_evals = 50), # specify 50 iterations
  search_space = search_space, # predefined hyperparameter search space
  tuner = mlr3tuning::tnr("random_search") # specify random search
)

AutoTuner -オブジェクトの train() -メソッドを呼び出すと、最終的にハイパーパラメータのチューニングが実行され、指定したパラメータに対して最適なハイパーパラメータの組み合わせが見つかる。

# hyperparameter tuning
set.seed(0412022)
autotuner_rf$train(task)
autotuner_rf$tuning_result
#>    mtry sample.fraction min.node.size learner_param_vals  x_domain regr.rmse
#> 1:    4             0.9             7          <list[4]> <list[3]>     0.375

15.4.2 予測マッピング

調整されたハイパーパラメータは、これで予測に使用することができる。 そのためには、フィットした AutoTuner オブジェクトの predict メソッドを実行するだけでよいのである。

# 最適なハイパーパラメータの組み合わせで予測
autotuner_rf$predict(task)
#> Warning: Detected version mismatch: Learner 'regr.ranger.tuned' has been trained
#> with mlr3 version '0.13.3', not matching currently installed version '0.14.0'
#> Warning: Detected version mismatch: Learner 'regr.ranger' has been trained with
#> mlr3 version '0.13.3', not matching currently installed version '0.14.0'
#> <PredictionRegr> for 84 observations:
#>     row_ids  truth response
#>           1 -1.084   -1.073
#>           2 -0.975   -1.050
#>           3 -0.912   -1.012
#> ---                        
#>          82  0.814    0.646
#>          83  0.814    0.790
#>          84  0.808    0.845

predict メソッドは、モデリングに使用されるすべてのオブザベーションにモデルを適用する。 モデリングで使用された予測因子として名付けられたラスタを含む多層 SpatRasterterra::predict() は空間予測、すなわち新しいデータに対する予測も行う。

pred = terra::predict(ep, model = autotuner_rf, fun = predict)
# hs = terra::shade(terra::terrain(dem, v = "slope", unit = "radians"),
#                   terra::terrain(dem, v = "aspect", unit = "radians"),
#                   10, 200) |>
#   terra::mask(terra::vect(study_area))
# pred = terra::mask(pred, terra::vect(study_area)) |>
#   terra::trim()
# library(tmap)
# pal = rev(hcl.colors(n = 15, "RdYlBu"))
# tm = tm_shape(hs) +
#   tm_grid(n.x = 3, n.y = 3, lines = FALSE) +
#   tm_raster(style = "cont", palette = rev(hcl.colors(99, "Grays")),
#             legend.show = FALSE) +
#   tm_shape(pred, is.master = TRUE) +
#   tm_raster(style = "cont", title = "NMDS1", alpha = 0.8,
#             legend.reverse = TRUE, palette = pal, midpoint    = NA) +
#   tm_shape(study_area) +
#   tm_borders() +
#   tm_layout(inner.margins = 0.02, legend.outside = TRUE)
# tmap_save(tm, "figures/15_rf_pred.png",
#           width = 12, height = 7, units = "cm")
knitr::include_graphics("figures/15_rf_pred.png")
植物相の勾配を予測するマッピングにより、はっきりとした植生帯が明らかになった。

FIGURE 15.5: 植物相の勾配を予測するマッピングにより、はっきりとした植生帯が明らかになった。

terra::predict() がモデル・アルゴリズムに対応していない場合でも、手動で予測を行うことができる。

newdata = as.data.frame(as.matrix(ep))
colSums(is.na(newdata))  # 0 NAs
# 0 があると仮定すると、より一般的なアプローチになる
ind = rowSums(is.na(newdata)) == 0
tmp = autotuner_rf$predict_newdata(newdata = newdata[ind, ], task = task)
newdata[ind, "pred"] = data.table::as.data.table(tmp)[["response"]]
pred_2 = ep$dem
# ここで、ラスタを予測値で埋める
pred_2[] = newdata$pred
# terra と我々の手動予測が同じかどうかをチェックする。
all(values(pred - pred_2) == 0)

予測地図には、はっきりとした植生帯( Figure 15.5)が描かれている。 ロマス山の植生帯の詳細については、Muenchow, Hauenstein, et al. (2013) を参照。 青い色調は、いわゆる Tillandsia-帯を表している。 Tillandsia は、特にロマス山脈の砂地やかなり砂漠的な麓で大量に見られる、高度に適応した属植物である。 黄色は草本植生帯で、Tillandsia-植生帯に比べ植物被度が高いことを表している。 オレンジ色はブロメリア帯を表し、種の豊富さと植物被覆率が最も高いことを特徴としている。 霧による湿度が最も高い気温逆転地帯(標高約750-850m)の真下で見られる。 逆転温度以上になると当然水分は減少し、再び砂漠化し、数種の多肉植物が見られるようになる(多肉植物帯;赤色)。 興味深いのは、空間予測によってブロメリア帯が途切れていることが明らかになったことである。これは、予測マッピングなしでは発見できなかった非常に興味深い発見であった。

15.5 結論

本章では、NMDS(Section 15.3)を用いて、ロマス モンゴル山の群集行列を順序付けした。 最初の軸は、調査地域の主な植物相の勾配を表し、部分的に R-GIS ブリッジ ( Section 15.2 ) を使って導き出した環境予測因子の関数としてモデル化された。 mlr3 パッケージは、ハイパーパラメータ mtry , sample.fraction および min.node.size ( Section 15.4.1 ) を空間的に調整するためのビルディングブロックを提供した。 調整されたハイパーパラメータ は最終モデルの入力となり、このモデルを環境予測変数に適用して植物相の勾配を空間的に表現した ( Section 15.4.2 ) 。 その結果、砂漠の真ん中にある驚異的な生物多様性を空間的に示すことができたのである。 ロマス山は絶滅の危機に瀕しているため、予測地図は保護区域を定める際の判断材料となり、地域住民に身近にあるユニークな存在であることを認識させることができるのである。

方法論の面では、いくつかの追加的な指摘ができる。

  • 2つ目の軸( )もモデル化し、2つの軸のモデル化されたスコアを1つの予測地図に統合して可視化する革新的な方法を見つけるのは興味深いことである。
  • もし、生態学的に意味のある方法でモデルを解釈することに興味があれば、おそらく(セミ)パラメトリックモデルを使うべきだろう (Muenchow, Bräuning, et al. 2013; A. Zuur et al. 2009; A. F. Zuur et al. 2017) しかし、少なくともランダムフォレストのような機械学習モデルの解釈を助けるアプローチは存在する(例えば、 https://mlr-org.github.io/interpretable-machine-learning-iml-and-mlr/ を参照)
  • 本章で使用したハイパーパラメータ のランダム化最適化よりも、逐次モデルベース最適化(SMBO)の方が望ましいかもしれない (Probst, Wright, and Boulesteix 2018)

最後に、ランダムフォレスト と他の機械学習 モデルは、多くのオブザベーションと多くの予測因子、この章で使われるよりもはるかに多く、どの変数と変数の相互作用が応答を説明するのに寄与するかが不明である設定で頻繁に使用されることに注意しておこう。 さらに、その関係は高度に非線形である可能性もある。 私たちのユースケースでは、レスポンスと予測変数の関係はかなり明確で、非線型はわずかであり、オブザベーションと予測変数の数は少ない。 したがって、線形モデル を試してみる価値はあるかもしれない。 線形モデルは、ランダムフォレスト( )モデルよりも説明や理解がしやすいので好まれ(パーシモンの法則)、さらに計算負荷が少ない(演習を参照)。 線形モデルがデータに存在する非線形性の程度に対処できない場合、一般化加法モデル (GAM) を試してみることもできる。 ここで重要なのは、データサイエンティストのツールボックスは複数のツールで構成されており、目の前のタスクや目的に最適なツールを選択するのはあなたの責任であるということである。 ここでは、ランダムフォレストのモデリングと、それに対応した空間予測への利用方法を紹介したいと思われる。 この目的のためには、反応と予測因子の関係が既知の、よく研究されたデータセットが適切である。 しかし、これはランダムフォレストモデルが予測性能の面で最良の結果を返したことを意味するものではない(演習を参照)。

15.6 演習

The solutions assume the following packages are attached (other packages will be attached when needed):

E1. Run a NMDS using the percentage data of the community matrix. Report the stress value and compare it to the stress value as retrieved from the NMDS using presence-absence data. What might explain the observed difference?

E2. Compute all the predictor rasters we have used in the chapter (catchment slope, catchment area), and put them into a SpatRaster-object. Add dem and ndvi to it. Next, compute profile and tangential curvature and add them as additional predictor rasters (hint: grass7:r.slope.aspect). Finally, construct a response-predictor matrix. The scores of the first NMDS axis (which were the result when using the presence-absence community matrix) rotated in accordance with elevation represent the response variable, and should be joined to random_points (use an inner join). To complete the response-predictor matrix, extract the values of the environmental predictor raster object to random_points.

E3. Retrieve the bias-reduced RMSE of a random forest and a linear model using spatial cross-validation. The random forest modeling should include the estimation of optimal hyperparameter combinations (random search with 50 iterations) in an inner tuning loop (see Section 12.5.2). Parallelize the tuning level (see Section 12.5.2). Report the mean RMSE and use a boxplot to visualize all retrieved RMSEs. Please not that this exercise is best solved using the mlr3 functions benchmark_grid() and benchmark() (see https://mlr3book.mlr-org.com/perf-eval-cmp.html#benchmarking for more information).