R

フレームワーク 事前準備：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構文の中身は次のように展開できる（下付き文字がイテレータに対応している）。...

効果量としての標準化偏回帰係数 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....