RでGIS GIS用のRパッケージ sfによるベクター処理 データを読み込む フィールド編集 レイヤー間の紐づけ まとめ RでGIS GISというとArcGISもしくはQGISを想定しがちだけど、個人的にはこれらGUIに依拠したSoftware（マウスでカチカチする系）は好みではないし、おすすめしない。というのも、何度もマウスでクリックしながら作業を進めるので、以下のリスクがある。 作業行程の再現性を担保しにくい エラーがあった場合、すべての行程を繰り返す必要がある うっかりミスしやすい これらの理由から、私はGISもスクリプトベースで作業行程を記録すべきと思う。スクリプトベースのGISというとPostGISが標準だったと思うが、最近ではこれらの関数群はほぼすべてRで使えるようになっているし、工夫すれば計算速度もArcGISと遜色ない。しかし、スクリプトベースのGISに関する資料は英語によるものがほとんどのため、アクセスしにくいといえばそうかもしれない。そこで、RのGIS処理について、いくつかの記事にわけて紹介したいと思う。なお、このブログで書かれていることの大半はTaro Mienoさん(University of Nebraska Lincoln)のオンライン資料で勉強させていただいた。とても事細かな説明があるので、この記事を読んで興味がわいた方はぜひこちらで勉強するといいと思う。 GIS用のRパッケージ どのようなタイプのレイヤー（ベクター、ラスター）を扱うかによって、必要となるパッケージは変わる。基本的な作業はsf, raster, terraあたりで済むが、多少込み入った作業には他の補助的なパッケージも必要になる。以下の表に簡単にまとめる。 パッケージ名 レイヤタイプ 備考 sf ベクター ベクター処理の場合はsf一択。大体これで足りる（Website）。 rmapshaper ベクター フィーチャーが重いときに頂点を削って計算を高速化する。（Website） raster ラスター ラスター処理のデファクトスタンダードだったが、C++ベースのterraに移行中（Website）。 terra ラスター ラスター処理の新しいデファクトスタンダードになりつつある（Website）。 stars ベクター、ラスター 時空間情報のついたラスター処理に特化。ラスター｜ベクター変換の際にも重宝する（Website）。 exactextractr ベクター、ラスター フィーチャーごとにラスターデータを集計するときに使う（Website）。 tmap NA GISレイヤーの描画（Website） sfによるベクター処理 データを読み込む ベクター処理の鉄板であるsfを使ってみる。sfについてくるデータセットを使う。拡張子を見てもらうとわかるが、呼び出されるファイルはシェイプファイル（....

応答変数（Y軸）の異なる図を並べる 軸のスケールをパネルごとに変える 変数変換 特殊文字 軸名 軸の値 パネルストライプ 見栄えのいい図表を作ることはとても好きで、不必要なほどにこだわってしまうこともある。しかし、「査読コメントに対応するために、図表を作り直すこと」は大嫌いであった*。これには理由がある。学生のころはＲの図表作成能力が低かったため、Rで図表のベースをつくったら、細かい調整をパワポやイラストレーターでしていたのだ。この作業は馬鹿にならない時間がかかるのだが、ちょっとした解析の修正や、リバイスの度にやり直しになる。 この作業による時間のロスが無駄だと感じたため、そのまま出版できる図表をコードを走らせるだけで作れるようせっせと豆知識をためてきた。最近ではRmarkdownと合わせれば、MicrosoftのOfficeに頼らずともすべての作業がRで完結する。ここでは、ggplot関連で案外わかるまで時間のかかったものに焦点をあててまとめる。もっといい書き方もあるかもしれないので、そのときはこっそり教えてほしい。今回は種数と面積の関係を模した以下のダミーデータを使って例を示す。 # dummy data x <- runif(100, 0.1, 1000) # hypothetical area m <- model.matrix(~log(x)) # model matrix y <- rpois(length(x), exp(m %*% c(log(5), 0.5))) # hypothetical richness df0 <- tibble(area = x, gamma = y, group = rep(letters[1:4], each = 25)) %>% mutate(alpha = rbeta(length(y), 5, 5) * gamma, beta = gamma / alpha) %>% pivot_longer(cols = c(alpha, beta, gamma), names_to = "metric", values_to = "diversity") print(df0) ## # A tibble: 300 × 4 ## area group metric diversity ## <dbl> <chr> <chr> <dbl> ## 1 526....

Package tidyverse patchwork sf/raster/stars/exactextractr whitebox Shortcut Code block label (Ctrl + Shift + R) Multi-line (un)comment (Ctrl + Shift + C) Pipe (Ctrl + Shift + M) Package tidyverse データ整理や図表作成に有用なさまざまなパッケージをまとめたもの。もはやこのパッケージなしにはRを使えない。とくにdplyr およびggplot2 に含まれる関数群にはお世話になりっぱなしである。このあたりの関数の使い方はWeb上にあふれているので割愛。dplyrであればHeavyWatalさんのWebsiteが分かりやすい。ggplot2に関してはFrom data to Vizがビジュアルから入れるのでとっつきやすい。 patchwork ggplot2とセットで使うとことを想定したパッケージ（patchwork）。データフレームをグループごとにプロットする場合、facet_wrapやfacet_gridなどの関数を使うことが多いと思う。これらは非常に有用な関数なのだが、すべてのパネルで同じ構造（xとyが一緒など）をとる必要がある。しかし、論文の図を作る時、フォーマットの異なる図を横に並べて一つの図としてまとめたいことも多いと思う（例えば散布図と箱ひげ図を並べる、など）。そんなときに役立つのがpatchworkである。このパッケージを使うと、複数のggplotオブジェクトを好きなように配置できる。irisを使って例を下に示す。 pacman::p_load(tidyverse, patchwork) # scatter plot g1 <- iris %>% ggplot(aes(x = Sepal.Width, y = Sepal.Length, color = Species)) + geom_point(alpha = 0.3) + theme_bw() # box plot g2 <- iris %>% ggplot(aes(x = Species, y = Sepal....

Structure Introduction Methods Results Discussion Additional point Paragraph Sentence 論文のドラフトをあげたら、以下の点を確認する。 Structure Introduction パラグラフ間の論理的つながりが明確 大きな枠組みの中で新奇性を述べている 対象とするシステムの「必然性」が明確 データの質と整合性が取れている 仮説が単純かつ明快 Methods 情報を出す順番が適切 パラグラフごとに情報が整理されている 理論・統計は専門外の人にも伝わるように書かれている パラメータ名の重複や説明漏れがない パッケージやソフトウェアのVersion情報に誤りがない Results Methodsを詳しく読まずとも結果が理解できる 結果の意味するところが述べられている 不確実性の記述が適切 効果量の記述が明確 Discussion 仮説・問いに対する答えが明確であり、答えの意義を広い枠組みのなかで捉えている（最初のパラグラフ） 各結果の解釈が端的に述べられている 複数の解釈がある場合、それらが網羅されている 弱みや限界に言及がある 今後の大きな展開がイメージできる（最後のパラグラフ） Additional point 読み落とされたくない重要情報は表現を変えて複数個所に配置する（MethodsとDiscussionなど）...

一見、線のようにみえる川ですが、上から俯瞰すると様々な川が合流を繰り返し、まるで樹木のような形をした「水系網」を作り上げます。その成因について、地形学や水文学の観点から長く研究されてきました。一方、この水系網の樹状形は、そこにすむ生物の成り立ちとどんな関わりがあるのでしょうか？ オーストラリア中央部にみられる樹状の水系網。 Photo Credit: Google Earth https://www.google.co.jp/intl/ja/earth/ ここ数年、私はこのテーマにかかわる研究を進めて来ました。その結果、水系網の分岐構造が複雑になるほど、（１）そこにすむ生物の集団はより安定的に維持され（Terui et al. 2018, PNAS）、（２）流域全体の生物多様性は高くなり（Terui et al. 2021, PNAS）、（３）さらには食物連鎖の長短にまで波及する（Pomeranz et al. 2021, preprint）ことが分かってきました。複雑な水系網では環境が多様なため、さまざまな生物を支える基盤が整っている、というのがおおざっぱな説明です。手法としては、自然界で起きているであろう生物と環境の相互作用を数式として記述し、そこから導かれる予測を現地データで検証しています。 いまになって論文化していますが、これらの発想はいずれも、博士課程在学時（2011-2014年）に川を泳いでいるときに思いついたものです。私の博士課程の研究は、川にすむカワシンジュガイという生物の空間的な分布を調べ、その成因を探るものでした（専門的にはメタ個体群構造を調べる；Terui et al. 2011, 2014）。その過程で、北海道の朱太川を4年かけて泳ぎまわりました（合計90km）。それほどアウトドアな人間ではない私は、ひとつの流域の中にある川が、これほど多様であることを知りませんでした。泳ぎながら貝を探していく中で、いくつもの川との合流地点を通り過ぎるのですが、合流を境に急に水が冷たくなったり、逆に温かくなったり、ときには底石の感じが変わったりするのです。川によって水の由来（湧水かそうでないか、など）や石の供給源（山）が変わるからなのですが、この調査をするまで気にも留めていませんでした。このとき、多様な川が織りなす「水系網」に強い興味がわきました。水系網の複雑さは、生物多様性を支えるとても重要な要素なのでは、と（これは川をよく見ている方なら感覚的にわかっていることと思います）。 2011年の北海道・朱太川における調査。写真中央は共同研究者であり、魚類研究者の宮崎博士。 しかし、当時はこの発想を科学へと昇華させることができませんでした。なぜなら、数理モデルやビッグデータを扱う能力がなく、この感覚的な仮説を科学の言葉で表現できなかったからです。仮に私がそうした技術を持ち合わせていても、当時の知識の乏しさから生態学における適切な位置づけを与えることができず、有象無象の研究として埋もれていたことでしょう。このため、しばらくはこの発想もとに研究を進めることはできず、学位取得後のポスドク期の後半に至るまで（2014-2017年）、全く別の研究に時間を投資しました。 転機があったのは2017年です。ポスドクのプロジェクトで関わっていた研究者の方が、北海道の保護水面で20年近く魚類のモニタリングが行われていることを教えてくれたのです。このとき、長い間抱えていた妄想が現実味を帯びはじめました。このモニタリングでは、多数の流域で魚類群集の時間変化が調べられており、水系網の形の意味を調べる上でまさに理想的なデータだったからです。さらにポスドク中、以前から興味のあった数理モデルの勉強を独学で進めていたので、このころには自分でシミュレーションモデルを作れるようになっていました。水系網の研究分野では、生態学の理論的枠組みが整理されていなかったので、それならば「自分で数理モデルの枠組みを作り、そこから導かれる予測をこの現地データで検証しよう」と考えました。そこから数年の間に渡米したのですが、その後に出会ったアメリカの研究者の協力もあり、今では日米をまたがる研究へと進展しました。こうして振り返ってみると、「過去のとりとめのない妄想」と「後に全く別目的で得た研究技術」の邂逅が、最近の研究成果につながったと言えます。いつ、どこで何が役に立つのか、本当にわからないものです。（もちろん、研究の実現にあたって、私に足りない部分を補ってくれた共同研究者の方々の協力が欠かせなかったことは言うまでもありません。） 生態系の構造的複雑さに注目する研究者は世界でも数えるほどしかいませんが、今回の一連の研究が、この研究分野の発展につながればと思います。個人的には、全く別の生態系 ‐ 例えば樹木や草本の樹状構造 - でも、それらを生息場とする微生物などで似た現象が起きそう、と考えています。

パラメータセットを用意する foreachの導入 Plotを描く 前回はfunctionを使ってシミュレーションモデルを関数化する方法を書いた。今日はこの関数を使い、様々なパラメータの下でシミュレーションを効率的に走らせるコードを書いてみたい。 パラメータセットを用意する 以前作った個体群動態のシミュレーションモデルを、ここの例題として再度利用する： \[ ln~n_{t+1} = ln~\lambda + ln~n_t + \epsilon_t\\ \epsilon_t \sim N(0,\sigma_{\epsilon}^2) \] 前年の個体数\(n_{t}\)に集団増加率\(\lambda\)が掛け算され、翌年の個体数\(n_{t+1}\)が決まるが、そこにはランダムな環境変動の影響（\(\epsilon_t\)）もある、というものであった（上の式では対数スケールのため足し算になっている）。これをコードとして書き下し、関数化したものが以下： sim_geomodel <- function(n_step, lambda, sd_eps, n1 = 10) { log_n <- NULL log_n[1] <- log(n1) eps <- rnorm(n = n_step, mean = 0, sd = sd_eps) for(t in 1:(n_step - 1)) { log_n[t + 1] <- log(lambda) + log_n[t] + eps[t] } n <- exp(log_n) df_dynamics <- dplyr::tibble(n_step = 1:n_step, n = n) return(df_dynamics) } この関数sim_geomodelを使い、パラメータの値を変えた時に、50年後の個体数の予測値がどう変わるのか調べてみよう。パラメータを変えながらパターンを予測することで、どのパラメータが集団動態にどんな影響を及ぼすのかを調べることができる。今回の場合、パラメータは二つあるので（lambda, sd_eps）、これらのパラメータセットを作るところから始める。パラメータセットを作るには、パラメータ値のすべての組み合わせを考える必要がある。この場合、expand....

functionの導入 シミュレーションを関数化 前回はfor構文を使った至極簡単なシミュレーションモデルを作ってみた。しかし、中には「こんなめんどくさいスクリプトを毎回書くのか。。。？」などと思われた方もいると思う。そんなことはないので安心してほしい。一つのまとまった作業を関数化することで、スクリプトの量をかなり減らすことができる。 functionの導入 function関数を使うことで一度書いたモデルを使いまわすことができる。最初からシミュレーションモデルを関数化すると説明が煩雑になってしまうので、まずは変動係数CVを推定する関数mycv なるものを作ってみよう。まずは正規分布に従う乱数をrnormを使って生成する。 # 100 random values following a normal distribution set.seed(123) # for reproducibility y <- rnorm(n = 100, mean = 50, sd = 25) # show the first 10 elements print(y[1:10]) ## [1] 35.98811 44.24556 88.96771 51.76271 53.23219 92.87662 61.52291 18.37347 ## [9] 32.82868 38.85845 変動係数は以下のスクリプトで推定できる。 cv <- sd(y) / mean(y) print(cv) ## [1] 0.4366692 しかし、なんだか毎回二つの関数sd とmean を組み合わせて変動係数を計算するのは面倒くさい。なので、これらの作業を一挙にやってくれる関数をつくってみよう。 mycv <- function(x) { cv <- sd(x) / mean(x) return(cv) } x という引数に基づいて、SDを平均で割るという作業を自動的にやってくれる関数mycv を定義している。function() のカッコの中に引数として使いたい変数をいれておく。そうすると、そこに使ってほしい値をぶち込むと、関数内で定義された作業を自動的に行ってくれる。return のところでは、何を計算結果（返り値）として返してほしいかを指定している。早速mycv を使ってみる。以下では、x = y とし、関数内のx にy を「代入」している。...

なぜシミュレーション？ for loop 集団動態モデル ランダムネスを加える モデルの拡張 なぜシミュレーション？ 生態学に慣れ始めてきたころ、いわゆる「理論研究」と言われる類の論文も読み始めるようになった。最初は難解で何をしているのかわからなかったが、分かってくるととても力強いアプローチだなぁと感じるようになり、自分で作ってみたいと思うようになった。というのも、私はフィールドを中心に研究をしていたけれども、野外のデータはあまりにも雑多で、その解釈に困ることが多かったからだ。例えば、ある魚と別の魚が餌をめぐる競争関係に興味があり、「この二種は競争関係にあるので、一方の個体数が多い場所では、もう一方の個体数は少なくなる」という仮説を立てたとしよう。野外で両種の個体数の間に負の相関が認められたとしても、「おお仮説通りのパターンだ、競争に違いない！」と単純に喜ぶことはできない。同じパターンを生み出す仕組みがあまりにもたくさんあるからだ（両者の好きな環境が全く異なるだけかもしれない）。 こうした理由から、自分がフィールドで集めたデータをもとに論文を書くとき、（特にDiscussionで）もどかしい思いをする。思いっきり「これだ！」と断言したいのに、あれやこれやと言い訳しなければならないからだ。実験で検証可能な仮説ならば、実験するに越したことはない。しかし、見たい現象が生態系スケールとかになってくると、実験などほぼ不可能だ。できたとしても億単位の研究費が必要になる。 そんなとき、シミュレーションが役に立つ。ある仕組みをこちらで勝手に想定し、そこから導かれるパターンがどんなものかを見るのがシミュレーションだ（生態学の数理モデルが全般的にそうですが）。つまり、観察されたパターンから仕組みを推論する統計モデリングの全く逆のことをするといってもいい（Figure 1）。興味のある仕組み以外を排除あるいはコントロールできるので、その仕組みがどんな時にどんなパターンを生み出すのか知ることができる。 Figure 1: Conceptual diagram for the roles of theoretical and statistical models. Theoretical models (generally) predict patterns under certain mechanisms (and assumptions) while statistical models infer mechanisms behind observed patterns と、ここまでは論文を読んでいれば納得できるのだが、いかんせんどうやってスクリプトを書けばいいのかわからない…というのが学生のころの悩みだった。統計解析のリソースはオンラインにかなり落ちているので自分でいくらでも勉強できたが、シミュレーションモデルは本当にスクリプトのリソースが少ない。あったとしても、これからやろうとしている人向けには書かれていない。それが今回の（たぶんシリーズ的に）書こうと思っているポストのモチベーション。 for loop こまごましたことはあるのだが、まずはfor構文をつかって簡単なシミュレーションモデルを作ってしまおう。for構文とはなんぞや、という人もいるかもしれないので、ここで簡単に説明しておく。端的にいうと、コンピュータに「同じ作業を繰り返せ」と指令するコマンド。for(i in 1:3) { XXX }という形で書くのだが、これはiがイテレータと呼ばれるもので、繰り返しのユニットに対応するものだ。このスクリプトで言えば、「XXXという作業を1から3まで回してほしい」と指令している。これだとわかりにくいので、下の例をみてみよう。 # create a vector with 11, 12, 13 y <- NULL x <- 11:13 for(i in 1:3) { y[i] <- x[i] * 2 } ここではyとxというオブジェクトを作り、xのほうに11から13の値を代入している。やりたいことは、xの値をそれぞれ2倍することである。for構文の中身は次のように展開できる（下付き文字がイテレータに対応している）。...

効果量としての標準化偏回帰係数 Average predictive comparison Rで実装 効果量としての標準化偏回帰係数 線形回帰は、生態学の一般的な解析手法になっている。Rで多数の関数が用意されているため解析も容易で、回帰係数の95%信頼区間やp値などでその効果の有意性もカンタンに検討できる。一方、近頃では有意性だけでなく、「効果量（Effect size）」にも注目したほうがいいとの見方が広まっている。この効果量には様々な指標があるけれども、線形回帰の文脈でいえば標準化偏回帰係数がもっとも広く使われている。 そもそも偏回帰係数の意味とは何かというと、対応する説明変数xが単位量（つまり1.0）だけ増えたときに、応答変数がどれだけ変化するかを示している。しかし、説明変数間で単位*1が違うと（例えばある変数はm単位で計られているのに対し、別の変数は㎝単位で計られている）、偏回帰係数はその単位の違いの影響を直に受けてしまう。そこで、それぞれの説明変数をそのばらつき（標準偏差SD）で割り（標準化）、それら標準化された説明変数に対する偏回帰係数を推定する。標準化された説明変数の単位はそろっているので、推定された標準化偏回帰係数も比較できるものになっているはず、というものだ。 しかし、これには一つの問題がある。現代の一般的な統計モデルでは、効果量としての標準化偏回帰係数の解釈が「直感的ではない」のである。例としては以下のようなものがある。個体数のカウントデータを線形モデルで表現する場合、誤差構造としてはポアソン分布（あるいは負の二項分布）を用いることが多い。ポアソン分布の平均パラメータは負の値をとることができないので、Rではデフォルトで対数リンク関数が実装されている。この時、ポアソンモデルの偏回帰係数のなにがどう「直感的」でないのか、下の例を見ながら考えてみる。 まず、下記のスクリプトで簡単な仮想データを作る。 # simulated data set.seed(111) n_sample <- 100 x1 <- rnorm(n_sample, mean = 10, sd = 1) x2 <- rnorm(n_sample, mean = 10, sd = 3) X <- model.matrix(~ x1 + x2) b <- c(0.1, 0.2, 0.2) # true parameters y <- rpois(n_sample, exp(X %*% b)) # simulated y ポアソン分布を仮定したGLMを当てはめる。 \[\begin{equation} y_i \sim Pois(\lambda_i)\\ ln \lambda_i = \alpha + \beta_1 x_{1,i} + \beta_2 x_{2,i} \end{equation}\] # model fit; unstandardized m <- glm(y ~ x1 + x2, family = poisson) coef(m) ## (Intercept) x1 x2 ## 0....

理論と実証の統合 河川魚類群集のビッグデータ 生物の移動パターン +アルファ 連絡先 私の運営する研究室では、海外学振あるいは学振PD、DCの制度を利用した研究員の受け入れを積極的に行っていきたいと考えています。私の研究室では以下のような興味を持った方が特に恩恵を得られると思います。 理論と実証の統合 私の研究室では理論と実証の統合に重きを置きたいと思っています。以下のような興味をもつ人が特にマッチングがよいと考えています。（１）これまでフィールドばかりやってきたが、理論も取り入れてもっと説得力のある研究デザインを組みたい、（２）理論をやってきたが、これからは実証も見据えた理論研究を展開したい。 私がカバーしている理論研究の分野は、生態系におけるネットワーク構造と集団・群集構造の関係、食物連鎖長の制限要因、生物の移動とそれが支える個体群構造に関する理論、などですが、例えば進化的視点を取り入れられるような方が来てくれると私の視点も広がるのでとてもうれしいです。 河川魚類群集のビッグデータ 当研究室では、魚類群集の大規模データにアクセスがあります。例えば、日本の北海道全域およびアメリカ中西部において、20000点を超える場所で調査した魚類群集データ（*）を整備したものがありますので、このデータベースを利用した研究計画を立てることが可能です。時系列データに関しても、北海道であれば100地点以上における20年にわたる魚類群集の継続的モニタリングデータがありますし、昨年には共同研究を通じてグローバルスケールで時系列データを収集したデータベースも出版しています（RivFishTIME）。これらのデータのシェアを通じ、国内外の共同研究者とも積極的に議論をしていく良い機会にもなると思います。また、これらのデータを用い、理論予測の検証に用いることも可能だと思います。 ＊既存の文献や州政府からデータ提供を受けたものです（私もごく一部のデータは集めてますが）。貢献された方々に感謝。 生物の移動パターン データベースではなかなか得られない「生物の移動」に関するデータを現地で集めています。私の所属するノースカロライナ大学グリーンズボロ校には、サテライトキャンパス（車で20分ほど）に小河川が流れています。ここではちょうどいい種数（Sunfish 3種、Chub 2種が群集の8-9割を占める川）があり、これらの魚種の移動が何に影響されているのかを調べています。長期的なプランとしては、詳細な移動パターン（移動の個体差など）が集団動態の安定性あるいは多種共存にどう関わっているのかを追求したいと考えています。 このサイトを活かしたフィールド研究あるいは理論予測の実証などを特に歓迎しますが、ほかにもザリガニや巻貝などの生物も多数いますので、本人の興味に応じて研究プロジェクトを拡張することも可能です。そのほか、もう少し離れた場所であればもっと原生的な河川でも調査が可能ですので、時間に余裕をもって相談くだされば調整可能です。 +アルファ 補足的な要素として、私は海外学振からアメリカで就職した経歴がありますので、そのあたりの経験をシェアすることも可能です。文化的に日本とかなり異なるところが多いので、そういった意味ではサポートできるかもしれません。アメリカ国内の研究者とのコネクションづくりをお手伝いもできると思います。 連絡先 そのほか、河川-陸域の食物網のつながりや、水系網と鳥類の分布の関係などにも興味があります。もちろん、上記にないようなテーマでも、可能な限りサポートしたいと考えています。興味のある方は照井a_terui at uncg.eduまでご連絡いただけますと幸いです。