混合モデル(Mixing model)

地球上には様々な元素があり、その性質を決めるのは陽子数です。しかし、同じ元素でも(陽子数が同じ)、中性子数が異なるために性質は同じくしながら質量の異なる「安定同位体」が存在します。混合モデルは、安定同位体の比が捕食を通じて濃縮される特性を利用し、捕食者に対する餌生物の相対的貢献度を推定する手法です。炭素(C)や窒素(N)の安定同位体比がもっとも一般的で、それぞれの濃縮係数(栄養段階が一つ上がる際の安定同位体比の濃縮の程度)についての情報が充実しています。

古典的な混合モデル(IsoSourceなど)では平均値しか扱われてこなかったため、種内の安定同位体比のばらつきや濃縮係数の不確実性は無視されてきました。しかし最近では、これらの不確実性を明示的に取り入れたベイズ混合モデルが一般的になっています(SIARなど)。SIARに代表されるベイズ推定を利用した混合モデルでは、以下の点が古典的な混合モデルと異なります。

  • 種内(あるいは任意に設定した栄養ノード内)の安定同位体比のばらつきを考慮
  • 濃縮係数を確率変数とすることで、安定同位体比の濃縮プロセスの不確実性を考慮

SIARはRで実装されており、だれでも簡単に利用できるようになっています。実装の手順は以下の通りです(詳しい解説はSIARのマニュアルを参照)。

#1 パッケージインストール
install.packages("siar")
library(siar)

#2 消費者のデータ
>consumer #should be matrix
      d15N      d13C
1 2.765054 -24.72063
2 9.857108 -24.22871
3 3.562164 -25.20802
4 7.623977 -23.86547
5 2.181553 -22.99689

#3 エサ生物のデータ
>prey
    Source mean_d15N  sd_d15N mean_d13C  sd_d13C
1 detritus      2.23 2.565637    -25.22 0.932000
2 plankton      0.13 3.711336    -22.33 3.461702
3   shrimp      4.43 3.930261    -23.34 2.922652

#4 濃縮係数のデータ
>tef
    Source mean_d15N sd_d15N mean_d13C sd_d13C
1 detritus      3.54    0.74      1.63    0.63
2 plankton      3.54    0.74      1.63    0.63
3   shrimp      3.54    0.74      1.63    0.63

#5 混合モデルの実行
result <- siarmcmcdirichletv4(consumer,prey,tef,concdep=0, iterations = 500000, burnin = 50000)

#6 推定結果の表示
Summary information for the output file ...
         Low 95% hdr High 95% hdr       mode      mean
detritus   0.1470012    0.9550798 0.43382742 0.5338503
plankton   0.0000000    0.5090717 0.03679839 0.2203550
shrimp     0.0000000    0.5362734 0.04320550 0.2457947
SD1        0.0000000    7.2432659 2.25476445 3.1661398
SD2        0.0000000    4.6901279 1.14716487 1.8523370

おまけ:グラフの表示

 siarproportionbygroupplot(model1) 

SPONSOR LINK

階層分割(Hierarchical variation partitioning)

階層分割(Hierarchical variation partitioning)は、線形モデルをもとにした分散分割の手法です。ある説明変数セットのもとで、すべての組み合わせのモデルを構築し、それらのフィッティング(決定係数など)を比較することで各説明変数の相対的な説明力を導きます。多重共線性の影響を受けにくいことから、野外データから各説明変数の相対的重要性を示すのに適した手法といえます。ただし、あまりにも説明変数が多い場合、モデル間の適合度指標の差を計算するプロセスが複雑化するため(ただ引いてるだけですが)、あまり信頼のおける結果は得られないように思います(説明変数はせいぜい3~5程度が限界?)。

階層分割はパッケージhier.partを用いることで容易に実装できます。

  • hier.part()
    パッケージ:hier.part
    family指定: 正規分布、ポアソン分布、二項分布
    ランダム効果:不可
    適合度指標:対数尤度(logLik)、決定係数(Rsqu)、Root-mean-square ‘prediction’ error(RMSPE)
    備考:Deviance explainedなどの適合度指標も計算は可能。
#sample
library(hier.part)
hier.part(y, X, family="gaussian", gof = "Rsqu")
# y: response variable
# X: dataframe of explanatory variables

# 結果の表示---------------------------------------------
# GOF for each candidate model
$gfs
[1] 0.00000000 0.31519073 0.01143161 0.32120974

# Independent and joint contributions of each variables
$IJ
             I         J      Total
x1 0.312484426 0.0027063 0.31519073
x2 0.008725311 0.0027063 0.01143161

# Percentage of independent contributions
# Note: the denominator is "total variance explained"
$I.perc
           I
x1 97.283609
x2  2.716391

SPONSOR LINK

線形モデル

glm()

パッケージ:デフォルト
family指定:正規分布、ポアソン分布、二項分布
ランダム効果:可
モデル選択:dredge, stepAIC
備考:特になし

# sample
glm(y ~ x, data, family = poisson)

glmer()

パッケージ:lme4
family指定:ポアソン分布、二項分布(正規分布の場合はlmer();R ver 3.Xあたりから仕様が変更?)
ランダム効果:可
モデル選択:dredgestepAIC
備考:
ランダム効果が多くなると推定精度が低下することが指摘されている。Bolkerさんの論文によれば、ランダム効果が3つ以上の場合はMCMCなど別の手法を採用したほうがいいとしている。ただし、これまで実装されていたmcmcsamp()という関数に関しては、分散の推定値がゼロ付近のときに適切な推定を行ってくれない可能性があるらしく、2017年5月現在はパッケージから関数が除外されている模様。なお、推定が収束しなかった場合、control=glmerControl(optCtrl=list(maxfun=20000))などと付け加えることで、試行数の上限を増やすことができる。

# sample 1: varying intercept model
glmer(y ~ x + (1|group), data, family = poisson)

# sample 2: varying intercept and slope model
glmer(y ~ x1 + (1 + x1|group), data, family = poisson)

# sample 3: nested random effect model
glmer(y ~ x1 + (1|region/group), data, family = poisson)
# region is a larger set of groups

# sample 3: nested random effect model
lmer(y ~ x1 + (1|region), data, control=lmerControl(optCtrl=list(maxfun=20000)))
glmer(y ~ x1 + (1|region), data, family=poisson, control=glmerControl(optCtrl=list(maxfun=20000)))
# when did not reach convergence with default settings

glm.nb()

パッケージ:MASS
family:負の二項分布に固定
ランダム効果:不可
モデル選択:dredge, stepAIC
備考:特になし

# sample
glm.nb(y ~ x, data)

glmmadmb()

パッケージ:glmmADMB
family:ポアソン分布、負の二項分布、ベータ分布、ベータ二項分布
ランダム効果:可
モデル選択:不可(関数を自作すれば可)
備考:zeroInflation=TRUEとすることでゼロ過剰モデルも実装可

# sample 1: poisson model with a random effect
glmmadmb(y ~ x + (1|group), data, family = "poisson")

# sample 2: negative binomial model with a random effect
glmmadmb(y ~ x + (1|group), data, family = "nbinom")

# sample 3: beta model with a random effect
glmmadmb(y ~ x + (1|group), data, family = "beta")

# sample 4: zero inflated poisson model with a random effect
glmmadmb(y ~ x + (1|group), data, family = "poisson", zeroInflation=TRUE)