HSM9_WM
WataruMKIMN
“2021-02-25”
Abstract
第Ⅲ部の中でも簡単な手法であるエンベロープアプローチと距離ベースアプローチの紹介第Ⅲ部イントロ
第Ⅲ部は実際に正則的ちを予測するために用いられる様々なアプローチを紹介している。 生態学的モデルは以下の3つに大別できる。 1. 記述的アプローチ：種の出現（在不在; 分布など）がどのような環境要因によって決まっているかを探求する。 2. 説明的アプローチ（仮説演繹的アプローチ）：ある生物学的知見を仮説として設定し、その仮説が支持/棄却できるかを確かめる。 3. 予測的アプローチ：予測力を最大限に高めるようなモデルを作成することを目的とする。（イメージとしてはDL研究などが近いかもしれない。DLと同様に中身がブラックボックス化してしまう恐れなどがある。）
説明力と予測力の間にはトレードオフがあることが知られている。
第Ⅲ部で説明する統計アルゴリズムはそれぞれ長所と短所があるため、モデル選択やマルチモデル推論などを行うと良い。
- マルチモデル推論：複数の統計モデルから一般化を行う作業のこと (Garamszegi & Mundry, 2014)。 “… it is best to employ an IT model averaging approach, a procedure that accounts for model selection uncertainty to obtain robust parameter estimates or predictions.” - Grueber et al, 2011
9章
在のみデータを用いた解析手法は非常に単純な規則に従っている。また偽不在データやバックグランドデータを作成する必要もない。 大別すると - エンベロープアプローチ - 距離ベースアプローチ
となる。
9.2 エンベロープアプローチ
地理座標エンベロープと環境エンベロープに分けられる、
- 地理座標エンベロープ：種の地理座標を全て含む最小単位の範囲を種の分布域と定義する。
- 環境エンベロープ：地理座標エンベロープに比べると、種の分布に影響を及ぼす潜在的な環境要因を考慮している。BIOCLIMなど。
エンベロープの欠点として、 - 種の分布の密度などを考慮しない（つまり、分布域の端と分布の中央が同じ意味を持つことになる！） - 環境要因に重み付けができない（つまり、全ての環境がエンベロープに対して同じ効果量を持つと仮定している！）
library(tidyverse)## ─ Attaching packages ──────────────────── tidyverse 1.3.0 ─## ✔ ggplot2 3.3.3 ✔ purrr 0.3.4
## ✔ tibble 3.0.6 ✔ dplyr 1.0.4
## ✔ tidyr 1.1.2 ✔ stringr 1.4.0
## ✔ readr 1.4.0 ✔ forcats 0.5.1## ─ Conflicts ───────────────────── tidyverse_conflicts() ─
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()library(biomod2) # install.packages("biomod2", repos="http://R-Forge.R-project.org")## biomod2 3.4.6 loaded.
##
## Type browseVignettes(package='biomod2') to access directly biomod2 vignettes.## Load the species and environmental datasets
mammals_data <- read_csv("hsdm/data/tabular/species/mammals_and_bioclim_table.csv") %>% dplyr::select(-1)## Warning: Missing column names filled in: 'X1' [1]##
## ─ Column specification ────────────────────────────
## cols(
## X1 = col_double(),
## X_WGS84 = col_double(),
## Y_WGS84 = col_double(),
## ConnochaetesGnou = col_double(),
## GuloGulo = col_double(),
## PantheraOnca = col_double(),
## PteropusGiganteus = col_double(),
## TenrecEcaudatus = col_double(),
## VulpesVulpes = col_double(),
## bio3 = col_double(),
## bio4 = col_double(),
## bio7 = col_double(),
## bio11 = col_double(),
## bio12 = col_double()
## )mammals_data## # A tibble: 8,542 x 13
## X_WGS84 Y_WGS84 ConnochaetesGnou GuloGulo PantheraOnca PteropusGigante…
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 -86.2 82.8 0 0 0 0
## 2 -84.8 82.8 0 1 0 0
## 3 -83.2 82.8 0 1 0 0
## 4 -81.8 82.8 0 1 0 0
## 5 -80.2 82.8 0 1 0 0
## 6 -78.8 82.8 0 1 0 0
## 7 -77.2 82.8 0 1 0 0
## 8 -75.8 82.8 0 1 0 0
## 9 -74.2 82.8 0 1 0 0
## 10 -72.8 82.8 0 1 0 0
## # … with 8,532 more rows, and 7 more variables: TenrecEcaudatus <dbl>,
## # VulpesVulpes <dbl>, bio3 <dbl>, bio4 <dbl>, bio7 <dbl>, bio11 <dbl>,
## # bio12 <dbl>この哺乳類データの中からアカギツネのデータを用いて種分布域エンベロープ (Species Range Envelope: SRE)を行う。
pred_BIOCLIM = sre(Response = mammals_data$VulpesVulpes,
Explanatory = mammals_data[,c("bio3", "bio7", "bio11", "bio12")],
NewData = mammals_data[,c("bio3", "bio7", "bio11", "bio12")],
Quant = 0)
pred_BIOCLIM_025 = sre(Response = mammals_data$VulpesVulpes,
Explanatory = mammals_data[,c("bio3", "bio7", "bio11", "bio12")],
NewData = mammals_data[,c("bio3", "bio7", "bio11", "bio12")],
Quant = 0.025)
pred_BIOCLIM_05 = sre(Response = mammals_data$VulpesVulpes,
Explanatory = mammals_data[,c("bio3", "bio7", "bio11", "bio12")],
NewData = mammals_data[,c("bio3", "bio7", "bio11", "bio12")],
Quant = 0.05)
par(mfrow=c(2,2))
level.plot(mammals_data$VulpesVulpes,
XY=mammals_data[,c("X_WGS84", "Y_WGS84")] %>% as.matrix(),
color.gradient = "grey",
cex=0.3,show.scale=F,
title="Original data")
level.plot(pred_BIOCLIM,
XY=mammals_data[,c("X_WGS84", "Y_WGS84")] %>% as.matrix(),
color.gradient = "grey",
cex=0.3,show.scale=F,
title="BIOCLIM 100%")
level.plot(pred_BIOCLIM_025,
XY=mammals_data[,c("X_WGS84", "Y_WGS84")] %>% as.matrix(),
color.gradient = "grey",
cex=0.3,show.scale=F,
title="BIOCLIM 97.5%")
level.plot(pred_BIOCLIM_05,
XY=mammals_data[,c("X_WGS84", "Y_WGS84")] %>% as.matrix(),
color.gradient = "grey",
cex=0.3,show.scale=F,
title="BIOCLIM 95%")
par(mfrow=c(1,1))左上から右下に向かって - 元データ - データの100%を用いたSRE - データの97.5%を用いたSRE - データの95%を用いたSRE
となっている。97.5%や95%などはエンベロープの感度を落としていることを意味する。
BIOCLIMは非常に単純であるが故に予測精度はよくない。しかし、単純なので基準予測として用いられることができる。
9.3 距離ベースアプローチ
単純なエンベロープアプローチの改良版という立ち位置。分析範囲全体を考慮する点が異なり、環境要因の重み付けなどに対応している。 いくつかの手法が開発されている。
- 主成分分析（PCA-sp）
- マハラノビス距離（MD）
- 生態的ニッチ要因解析（ENFA）
ここではENFAについて説明する。
ENFA
ENFAの大まかな流れは 1. 分析範囲全体の多次元環境の平均と、分析範囲（種の在する空間）の多次元環境の平均の差をとる 2. 分析範囲全体の多次元環境の分散と、分析範囲（種の在する空間）の多次元環境の分散の差をとる 3. 1.2.の差から調査地全体の環境分布とその種の環境空間上の分布を比較し、その種の生態ニッチを推定する
である (参考資料)。また、PCA同様に環境要因を無相関変数したりする。 以下でENFAを動かしてみる。
library(adehabitatHS)## Loading required package: sp## Loading required package: ade4## Loading required package: adehabitatMA## Registered S3 methods overwritten by 'adehabitatMA':
## method from
## print.SpatialPixelsDataFrame sp
## print.SpatialPixels sp## Loading required package: adehabitatHR## Loading required package: deldir## deldir 0.2-9 Nickname: "Morpheus and Euripides"##
## Note 1: As of version 0.2-1, error handling in this
## package was amended to conform to the usual R protocol.
## The deldir() function now actually throws an error
## when one occurs, rather than displaying an error number
## and returning a NULL.
##
## Note 2: As of version 0.1-29 the arguments "col"
## and "lty" of plot.deldir() had their names changed to
## "cmpnt_col" and "cmpnt_lty" respectively basically
## to allow "col" and and "lty" to be passed as "..."
## arguments.
##
## Note 3: As of version 0.1-29 the "plotit" argument
## of deldir() was changed to (simply) "plot".
##
## See the help for deldir() and plot.deldir().## Loading required package: adehabitatLT## Loading required package: CircStats## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select## Loading required package: boot##
## Attaching package: 'adehabitatLT'## The following object is masked from 'package:dplyr':
##
## idlibrary(pROC)## Type 'citation("pROC")' for a citation.##
## Attaching package: 'pROC'## The following objects are masked from 'package:stats':
##
## cov, smooth, varpc <- dudi.pca(mammals_data[,c("bio3", "bio7", "bio11", "bio12")], scannf=FALSE, nf = 2)
en <- enfa(pc, mammals_data$VulpesVulpes, scan=FALSE)
par(mfrow=c(2,2))
barplot(en$s) # the specialization diagram
scatterniche(en$li,mammals_data$VulpesVulpes, pts=T) # plot the niche
s.arrow(cor(pc$tab,en$li)) # meaning of the axes
左上からZ字に 1. 特殊化（specialization）のeigen value （上記資料ではS）：今回はbio3-bio12の4変数なので3本（次元圧縮はPCAと同様。） 2. specialization 1とmarginalityの散布図：環境全体が灰色で種が在する範囲は黒。 3. bio3-bio12のPCAで調整された値(?)とENFAで求めたspecialization 1とmarginalityのcorrelation
ENFAとBIOCLIMを比較してみようとすると、どちらかの値に範囲を揃える必要がある。 今回はENFAの値を在不在の2値情報に変換して比較する。変換のために{pROC}のroc()を用いた。
par(mfrow=c(2,2))
level.plot(mammals_data$VulpesVulpes, XY=mammals_data[,c("X_WGS84", "Y_WGS84")] %>% as.matrix(),
color.gradient = "grey", cex=0.3,show.scale=F, title="Original data")
level.plot(en$li[,1], XY=mammals_data[,c("X_WGS84", "Y_WGS84")] %>% as.matrix(),
color.gradient = "grey", cex=0.3,show.scale=F, title="ENFA")
roc_enfa <- roc(mammals_data$VulpesVulpes, en$li[,1])## Setting levels: control = 0, case = 1## Setting direction: controls < casesthreshold_enfa <- pROC::coords(roc_enfa, "best", ret=c("threshold"))
# Pred01=as.numeric(en$li[,1] > threshold_enfa)
Pred01 <- en$li[,1] %>% as_tibble() %>% mutate(value2=if_else(value>threshold_enfa %>% pull(threshold), value, 0)) %>% pull(value2)
level.plot(Pred01, XY=mammals_data[,c("X_WGS84", "Y_WGS84")] %>% as.matrix(), color.gradient = "grey", cex=0.3,show.scale=F, title="ENFA binary")
なんか3枚目（左下）が教科書と異なるような気がする。roc()部分とPred01=as.numeric(en$li[,1] > threshold_enfa)がわかっていない。