シャコ・エビ日記

シャコパンチ、エビパチン研究者の記録

カイカムリのスポンジ選択で学ぶ階層ベイズモデリング

このプレプリントで、カイカムリにサイズの異なる三つのスポンジを与えて、どのサイズを気に入るかテストしています。その解析部分をここで解説してみます。それぞれの選択肢を選ぶ確率が身体の大きさと脚の欠損度合いとどのような関係にあるか、統計モデルをつくって推定してみようと思います。また、この選択行動に個体差があると考えてそれがどれくらいか推定してみます。松浦氏のStanとRでベイズ統計モデリング階層ベイズモデルとWAICに多くを負っています。

https://camo.githubusercontent.com/36952df365d9e5d79b6f29fd9b5a9eed2722332d/68747470733a2f2f64326d787565667165616137736a2e636c6f756466726f6e742e6e65742f735f324134354241463833433634453335323233303741383235304245423630454544384139463232343535443545424441463741433642303036413043363442315f313532363731353630303439305f4669675f6d6174657269616c735f5f6d6574686f64732e706e67

モデル構成

Sサイズを選んだ個体はいなかったので、MとLサイズを選んだ場合と、いずれも選ばなかった場合を考えます。Mの選択を基準とした、ソフトマックス回帰を行います。また、Mと無選択の両方を示した個体は1個体だけだったため、個体差はLの選択のみ考えます。甲幅Cwidth、脚の欠損度合いがLegLackです。

 \mu[n, 1] = 0
 \mu[n, 2] = bias\_l[ID[n]] + cwl * Cwidth[n] + ll * LegLack[n]
 \mu[n, 3] = bias\_no[n] + cwno * Cwidth[n] + lno * LegLack[n]
 n=1,..., N_{act}

 bias\_l[k] \sim Normal(bias\_l0, s\_bias\_l )
 k=1,...,N_{animal}

 Choice[n] \sim Categorical(softmax(\mu[n,])), n=1,..., N_{act}

Stanコード

gist.github.com

generated quantities ブロックで、functionsブロックで定義されている関数を実行し、モデルの予測力をWAICを用いて比較するための計算を行っています。ここでは、予測分布としては、新たな個体の新たな選択の予測を考えるのが適当ですから、個体差を表わすために導入されているローカルなパラメータ( bias\_l) を積分しています。高速化するために、階層ベイズモデルとWAICを参考にして、分けて計算をしています。

Stanを実行するRコード

gist.github.com

パラメータの事後分布

最近見つけた、{ggdistribute}でプロットしてみました。脚の欠損度合いは関係なく、甲幅は選択に影響がありそうだということが分かります。
f:id:katsumushi:20180604001602p:plain

また、個体差を表現する bias\_lの事後分布を全個体分プロットすると以下のようになります。
f:id:katsumushi:20180604001617p:plain

プロット

以下のプロット関数を実行すると、次の図*1が表示されます。身体が大きくなると、Lサイズを選ぶ確率が上がり、9 cmあたりを超えると、どれも選ばなくなっているのが分かります。

https://camo.githubusercontent.com/96328b57d2b3785820c95c39e7afb194e63ef102/68747470733a2f2f64326d787565667165616137736a2e636c6f756466726f6e742e6e65742f735f324134354241463833433634453335323233303741383235304245423630454544384139463232343535443545424441463741433642303036413043363442315f313532363237383035313637305f4669675f6265686176696f72616c5f63686f6963652e706e67


gist.github.com

パラメータの50パーセンタイルを用いてモデルに基づいてデータを発生させて*2、それを stat_density2d で重ねて、濃いところでデータ発生確率が高い様子を表現しました。my_categorical_logit_rng 関数は、次を rstan::expose_stan_functions 関数で読み込んでから実行してください。

https://gist.github.com/kagaya/c726f16d82b80026bb1924b408a72b5c


データ点が濃いところで出ているようなので、モデルの予測がうまく行っていると言ってよいかと思います。また、個体38は小さいにもかかわらず、Lサイズを好んで選ぶような個性が出ているのですが、階層モデルがこれには過剰に適合していないように見えます。

WAIC の計算

gist.github.com

で計算できます。ここではひとつだけなので、比較はしませんが、詳しくはプレプリントを御参照されるか、次回の記事で比較してみようと思います。

*1:正確にはその原図

*2:ここはしかしフルベイズ、つまり、せっかく求めたパラメータの事後分布をフルにつかって予測分布をつくったほうがよいので、そこを論文でも改訂するつもりです。