Rで混合モデル (R Advent Calendar 2013 16日目)

R

R Advent Calendar 2013の16日目です．

「Rと同位体分析: siarパッケージの紹介」の内容を文章にして，補足情報を載せました．


背景

生態学の分野では，ある生物がさまざまな食資源をそれぞれどのくらいの割合で摂取しているかを知ることが，重要な研究課題のひとつになっています．この研究課題に対して，ここ数十年のあいだに，同位体分析と混合モデル (mixing model) によるアプローチがよく行なわれるようになってきました．


同位体分析と混合モデル

同位体分析についてすごーく簡単に述べます．ある生物の体組織の同位体比を測れば，その生物の摂取してきた食物の内容がおおまかにわかります．同位体比というのは，質量数の異なる同一元素の存在比のことで，炭素の場合は，^{13}C/^{12}Cとか^{14}C/^{12}Cとかになります．たとえば生態学でもっともよく使われる炭素・窒素の同位体分析を利用することで，ある生物が，海・淡水・陸のどこから食資源を得ているか，肉食・草食のどちらに近いか，といったことがわかります．なお，元素Xの質量数nの同位体比は\delta^{n}Xというように表記し，単位は千分率 (‰)，標準物質からの差分として表されます．


具体例を出すと，以下のようになります．アマモとアオサのみを食べているガチョウがいたとして (図1)，それらの食資源とガチョウ体組織の炭素同位体比 (\delta^{13}C値) を測って連立方程式を解けば，食物全体に占めるアオサ，アマモの相対的な寄与率を求めることができます．

同位体比の測定結果 δ13C (‰)
ガチョウ (同位体分別補正済み) -13.2
アマモ -11.2
アオサ -14.1


図1 ガチョウとその食資源の炭素同位体比


δ13C{ガチョウ} = a δ13C{アマモ} + b δ13C{アオサ}
a + b = 1

同位体比測定結果をもとににこの連立方程式を解くと，以下のようになります．
アマモ寄与率 a = 0.31
アオサ寄与率 b = 0.69


用いる同位体の指標を2つに増やせば，3つの食資源の寄与率を計算できます．先ほどの計算に，新たに窒素同位体比 (\delta^{15}N) の結果と食資源の草を加えて，アオサ・アマモ・草のみを食べているガチョウを考えてみます (図2)．

同位体比の測定結果 δ13C (‰) δ15N (‰)
ガチョウ (同位体分別補正済み) -13.2 6.8
アマモ -11.2 6.5
アオサ -14.1 9.8
草 -30.9 4.4


図2 ガチョウとその食資源の炭素・窒素同位体比


δ13C{ガチョウ} = a δ13C{アマモ} + b δ13C{アオサ} + c δ13C{草}
δ15N{ガチョウ} = a δ15N{アマモ} + b δ15N{アオサ} + c δ15N{草}
a + b + c = 1

連立方程式を解くと，以下のようになります．
アマモ寄与率 a = 0.81
アオサ寄与率 b = 0.13
草寄与率 c = 0.06


混合モデルとは，こうした，食資源と消費者の同位体比を元に，それぞれの食資源の寄与率を計算するモデルのことです．どうしてこんな簡単な計算にモデルが必要かというと，以下3つの大きな問題があるためです．

  • 利用可能な同位体比の指標 (せいぜい2-4種類) に対して，摂取割合を推定したい食資源の数が多い場合，連立方程式を解くことができなくなる．
  • 生物現象のため，食資源や消費者体組織の値にばらつきがあり，平均値で代表させた値の計算では不正確なことがある．
  • 上の例では示していないが，食資源の元素濃度や，消費者体内での同位体の扱いの違いなど，さらに考慮しなければならない変数がある．


こうした問題を解決するため，モンテカルロ法による解析，総あたり計算，濃度依存混合モデルなど，さまざまな解決案が提示されてきました．このあたりの事情は，Phillips (2012) をご参照ください．


ベイズ推定を利用した混合モデル"siar"

そんななか，2010年に，ベイズ推定を利用した強力なモデルが，Parnell et al. (2010) によって提案されました．このモデルは，RパッケージとしてCRANに公開され，Google Scholarで調べてみると，論文の引用数はすでに320を超えています (2013年12月16日現在: 同位体生態学の論文では多いほうだと思います)．


このパッケージsiar (Stable Isotope Analysis in R) の紹介をするのが本エントリの目的です (やっとRの話が出てきました！)．ユーザインタフェースがつくりこまれており，ふだんあまりRに触れない方でも，けっこう簡単に使えてしまいます．


特徴は以下の通りです．

  • (用いる同位体比の指標数 + 1) 以上の食資源からの寄与率であっても計算できる．
  • ベイズ推定のフレームワークを採用しているので，値のばらつきをモデルに組み込み，各食資源の寄与率を事後分布として推定できる．事前分布も設定可能．
  • 食資源の元素濃度を組み込んだモデルや，消費者体内での同位体の扱いの変動などにも対応可能．


モデルの中身については，元の論文 (Parnell et al. 2010) をご覧ください．仕組みは非常に簡単なので，同位体生態学以外の分野にも応用が効きそうです．


デモ

以下に，様子を載せておきます．

# まずインストール
install.packages("siar")
# パッケージの読み込み
library(siar)
# メインの関数を開いてみましょう
siarmenu()

そうすると，以下のような選択式のメニューがあらわれ，あとはデータを読み込んで，数字を選んでいくだけで，ベイズ推定を組み込んだ混合モデルを利用したひと通りの解析ができます．

------------------------------- 
Welcome to Stable Isotope Analysis in R version 4.1 
Author: Andrew Parnell, University College Dublin
Please report bugs to: Andrew.Parnell@ucd.ie
------------------------------- 

Useful: Press 0 at a prompt to return to the main menu or Esc to exit. 
 
The available options are: 

 1: Load in some data (run this first)
 2: Run SIAR for a single group
 3: Run SIAR for multiple groups
 4: Plot single group proportions
 5: Matrix plot of proportions
 6: Plot of proportions by source
 7: Save parameter output to a file
 8: Summary information and convergence diagnostics
 9: Demo (for first time users)
10: Exit

Selection:

はじめは9のデモを選んでみると，一連の流れがよくわかります．上記に示したガチョウの食性データはこのデモから引用したものです．


関数群

自分好みに変数や設定を変えたりしたい場合は，ヘルプと，以下のガイドが参考になります．パッケージsiarの作者らが作成した文書で，環境中に読み込まれる関数などについて説明しています．これらの関数を使えば，いちいち，関数siarmenuによるユーザインタフェースを通して解析する必要がなくなり，通常ようなコマンドラインからの操作ができます．


SIAR V4 (Stable Isotope Analysis in R) An Ecologist's Guide (PDF)


たとえば，デモで見たのと同じ操作が以下のコマンドで可能です．

# データ読み込み
data(geese1demo)
data(sourcesdemo)
data(correctionsdemo)
data(concdepdemo)

# MCMCを走らせる
siarresult <- siarmcmcdirichletv4(
  data = geese1demo,
  sources = sourcesdemo,
  corrections = correctionsdemo,
  concdep = concdepdemo)

# 生データのプロット (図3)
siarplotdata(siarresult)

# 混合モデル解析結果のプロット (図4)
siarhistograms(siarresult)



図3 ガチョウとその食資源同位体比のプロット



図4 各食資源の寄与率推定結果


まとめ

混合モデルは，同位体生態学の研究以外にも応用可能です．「何か異なるものが混ざり合ってひとつのものができていて，元の混合比率を知りたい」という課題に対して，ここで紹介したパッケージsiarは強力な解法を提供してくれます．なにか機会があったら，ぜひ使ってみてください．


引用文献

いずれもオープンアクセスです．


Parnell AC, Inger R, Bearhop S, and Jackson AL. 2010. Source partitioning using stable isotopes: coping with too much variation. PLoS One 5:e9672.


Phillips DL. 2012. Converting isotope values to diet composition: the use of mixing models. J Mammal 93:342-352.