R

中心極限定理。統計を使う研究に身を置く人であれば、一度は聞いたことがあるであろう。 はて、その意味するところはなんだろうか？ 「サンプルサイズ（標本数）が増えれば、どんな確率変数も正規分布で近似できる」 たぶん、生態学の界隈では、このような理解がまことしやかに蔓延っている（昔私も信じていた）。その結果としてこういった記述をみることがある。 「応答変数は個体数であるが、サンプル数が多いことから正規分布で近似する」 結論から言うと、できません。さて、どこが間違ってるのだろうか。 わかりやすい説明を試みるが、自信はない。 Rで試す中心極限定理 統計の表現になれのない人のために、表現について少し解説を加えておく。ここでは、正規分布ではない確率分布として、ポアソン分布を例に挙げようと思う。確率変数1がある確率分布に従うとき、それはニョロニョロ（Tilde, ～）をつかって表現される。例えば、確率変数Xが平均２のポアソン分布に従う場合、次のように書く： \[ X \sim \text{Poisson}(2) \] 試しにRでXを生成してみる（rpois()はポアソン分布に従う乱数を生成）。 # produce 3 samples of X that follows a Poisson distribution with a mean 2 (X <- rpois(n = 3, lambda = 2)) ## [1] 1 2 2 この例では、３サンプル生成してみると、1, 2, 2 となった。ここで、最初の文言に戻ってみよう。 「サンプルサイズ（標本数）が増えれば、どんな確率変数も正規分布で近似できる」 この言葉を字面通りに受け取ると、サンプルをどんどん増やしていけばXはいずれ正規分布に近づくように解釈できる。もしそうならば、サンプル数を増やすほど、Xの頻度分布はきれいな左右対称の鐘状になるはずだ。試してみよう。まずは100サンプル。 # histogram with 100 samples X <- rpois(n = 100, lambda = 2) plot(table(X)) うーん、歪んでいる。１万ならどうだ。 # histogram with 10000 samples X <- rpois(n = 10000, lambda = 2) plot(table(X)) そろそろわかると思うが、Xのサンプルサイズを増やしたところでXはポアソン分布のままである。そう、中心極限定理は、そもそもこんなことは言っていないので当然である。...

RでGIS terraによるラスター処理 データの読み込み ベクターレイヤと組み合わせる crop() ：必要な部分のみ切り抜く exact_extract() ：ポリゴンごとに値を集計 最後に RでGIS 先日、RでのGIS作業にむけた簡単なポスト（RでGIS：ベクター）を書いた（最近だと思っていたが実に二年以上たっている事実に驚愕）。今回はそのフォローアップとして、ラスターデータの取り扱いについても簡単に紹介してみたい。今回はsf、terra、exactextractr、tmapの４つを使うので、先にこれらを読み込んでおく。 pacman::p_load(sf, terra, exactextractr, tmap, tidyverse) terraによるラスター処理 Rによるラスターデータの解析は長らくrasterパッケージが使われてきたようだが（私も初期はこちらを利用）、最近ではterraが主流になっており、ほかのGIS系パッケージとも互換性が確保されてきた様子（いくつか対応していないものもあるようだが、細かく把握していない）。terraのほうが圧倒的に計算速度が速いので、こちらの例を挙げる。 データの読み込み データの読み込みにはterra::rast()を使う。今回、ここで使うラスターファイルはCHELSAの気候データ（赤道付近で1km程度のメッシュ）。CHELSA_bio1_1981-2010_V.2.1.tif（1981-2010年の平均気温データ）のアメリカ、ノースカロライナ（NC）周辺を切り出し、.../static/gis/airtemp_nc.tifとして保存しておいた（Linkからダウンロード可能）。私はこのファイルパスを指定して読み込むが1、使用者は各自の保存先のパスを指定すれば読み込める。 # read data # here::here() is a function that returns an absolute path rs_air <- rast(here::here("static/gis/airtemp_nc.tif")) print(rs_air) ## class : SpatRaster ## dimensions : 325, 1064, 1 (nrow, ncol, nlyr) ## resolution : 0.008333333, 0.008333333 (x, y) ## extent : -84.32514, -75....

フレームワーク 事前準備：Github/Gitの環境構築 事前準備：Projectをつくる 事前準備：Projectの内部構造 編集する コミットする プッシュする まとめ 再現性のあるコーディングはとても重要だが、整理の仕方をまとめたものがほとんどない。個人の好みが大きく出るところだが、自然選択を経て私が到達したレポジトリの作り方をまとめてみる。前提条件としてR Studioを使っていることを想定するが、そうでなくとも基本構造は使える。 フレームワーク ローカル（自分のパソコン）とオンラインレポジトリ（Github）を連携する環境を作る。変更履歴を残しながらオンラインレポジトリに保存することで、デバッグ（エラーがないかチェックすること）・再現・シェアしやすくなる。大枠としては以下： レポジトリ（コードなど）を編集する 変更履歴をGitに記録する（コミット） 編集をGithubに反映する（プッシュ） １に戻る 事前準備：Github/Gitの環境構築 Happy Git and GitHub for the useRがかなり丁寧にまとめてくれているので、こちらをぜひ一読してほしい。以下は最低限の環境構築の手順を説明する。 Githubでアカウントをつくり、オンラインでのストレージ先を確保する。このアカウント名はレポジトリのリンクに必ず含まれるので、極端に長いものや、大文字を含むものは避けたい。このアカウントにローカルで作成したレポジトリを追加する。基本、Github上でレポジトリは編集せず、ローカルで編集したものをオンラインにプッシュする形で更新する。 Gitをインストールする。ローカルで機能するもので、編集履歴を自動管理するための装置（バージョンコントロールシステム）である。Gitでは、コミットすることで編集履歴を記録してくれる。変更箇所は自動認知してくれるため、この作業が非常に楽になる。基本デフォルトの設定のままインストールすれば問題ないが、念のためこちらを読んでおくといいかもしれない。 Git clientをインストールする。Gitはターミナルから直接操作もできるが、ややとっつきにくい。そこで、操作しやすいGUIを提供してくれるGithub Desktopをここでは使うことにする。ほかの選択肢についてはこちらを参照するとよい。 事前準備：Projectをつくる Gitの環境構築が済んだらProjectを作る（メニューバーのFile > New Project > New Directory > New Project ）。このとき、以下のウインドウでCreate a git repositoryにチェックし忘れないようにする。このProjectに論文一つの必要情報をすべて入れる（生データ、コード、論文原稿）。 R Projectの作成画面。“Create a git repository”にチェックし忘れないようにする。なお、Gitがインストールされていないとこのチェックマークは選択不可。 プロジェクトは、初期ファイルとして以下のファイルを含んでいる。 project_name.Rproj: プロジェクトのメタデータファイルのようなもの。プロジェクトを開きたいときはこれを開く。プロジェクトが開かれている場合は、R Studio右上にproject_nameがみえる。開いていない場合はProject: (None)となる。...

Publication list Coauthor list I have been spending hours to list my pubs for my CV or co-authors for a grant proposal. But now, we can automate this process with R package scholar (see R documentation for general guidance). With some stringr and dplyr functions, it’s pretty easy to export a table of pubs/coauthors via Rmarkdown (if you are unfamiliar with Rmarkdown, see Rmarkdown cookbook). A procedure would be:...

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....

パラメータセットを用意する 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構文の中身は次のように展開できる（下付き文字がイテレータに対応している）。...