Rでシミュレーション2 - function

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 を「代入」している。...

April 11, 2021 · Akira Terui

Rでシミュレーション1 - for loop

なぜシミュレーション? 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構文の中身は次のように展開できる(下付き文字がイテレータに対応している)。...

March 29, 2021 · Akira Terui

回帰モデルの効果量

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

March 28, 2021 · Akira Terui