BLOG

半透明のシンボル

半透明のプロットを描く場合、シンボルの中身の色を調整できるシンボルを指定します。ここでは、シンボルの引数pchは21とし(円の中身の色指定ができるシンボル)、bgで中身の色を指定します。pch=21の場合、colはシンボルの縁どりの色を指定します。

rgb()はred,green,blueの三原色の割合を指定することで色調整する関数です。最初の三つの引数がそれぞれRGBに対応しています。すべて0とすると黒、1とすると白になります。alphaが透明度を指定する引数で、0(完全な透明)から1(完全な不透明)の値をとります。

#points()内の引数の値を調整
plot(y~x, type="n", ann=F, axes=F)
box(bty="l")
points(x, y, cex=1.5, pch=21, col=NA, bg=rgb(0,0,0,alpha=0.3)
axis(1); axis(2, las=2)
mtext("Standardized lake area", side = 1, line = 3)
mtext("Fish density", side = 2, line = 3)

 

SPONSOR LINK

プロットの調整

下のplotを使い、見た目の調整方法についてまとめます。私個人の意見では、一度まっさらなプロット領域を用意し、好みの設定でデータや軸を追加していくのがのちのち微調整しやすいスクリプトになると思います。

調整前
#1. バックグラウンドとなるプロットを作る
plot(y~x, type="n", ann=F, axes=F)
#type="n"でプロット領域に何も書かない
#ann=F で両軸ラベルの消去
#axes=F で両軸の消去

#2. 外枠の追加
plot(y~x, type="n", ann=F, axes=F)
box(bty = "l")
#bty="l"でL字型の外枠

#3. データポイントの追加
plot(y~x, type="n", ann=F, axes=F)
box(bty = "l")
points(x, y, cex=1.5)
#cexでポイントサイズの調整:2とするとデフォルトの二倍の大きさ

#4. 軸の追加
plot(y~x, type="n", ann=F, axes=F)
box(bty = "l")
points(x, y, cex = 1.5)
axis(1); axis(2, las = 2)
#axisの第一引数は軸の位置を指定:1=下横軸、2=左縦軸、3=上横軸、4=右横軸
#las=2とすることで軸の値を横表記

#5. 軸ラベルの追加
plot(y~x, type="n", ann=F, axes=F)
box(bty = "l")
points(x, y, cex=1.5)
axis(1); axis(2, las = 2)
mtext("Standardized lake area", side = 1, line = 3)
mtext("Fish density", side = 2, line = 3)
#sideはラベル位置の指定:1=下、2=左、3=上、4=右
#lineはラベルの軸からの距離(行数)

※ステップ1~6まで継ぎ足しで書いてます

調整後

今回使ったサンプルデータの生成

N <- 50; a <- 3.4; b <- 0.2
x <- rnorm(N,0,1)
sigma <- 0.3

mu <- y <- NULL
for(i in 1:N){
  mu[i] <- a + b*x[i]
  y[i] <- rnorm(1,mu[i],sigma)
}

SPONSOR LINK

混合モデル(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) 

階層分割(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

線形モデル

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)

結果の表示

RからBUGSソフトウェアを動かすのスクリプトを入力すると、次のような計算プロセスが表示され、postに計算結果が格納されます。

Compiling rjags model...
Calling the simulation using the rjags method...
Adapting the model for 100 iterations...
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
Burning in the model for 500 iterations...
|**************************************************| 100%
Running the model for 1500 iterations...
|**************************************************| 100%
Simulation complete
Calculating summary statistics...
Calculating the Gelman-Rubin statistic for 2 variables....
Finished running the simulation

次のコードを入力し、結果を確認します。

summary(post)
        Lower95    Median   Upper95     psrf
alpha -0.02896274 0.1647402 0.3390248 1.001341
beta   0.63089605 0.7559828 0.9100826 1.001867

Medianが推定結果の中央値で、Lower95とUpper95がいわゆる95%信用区間(CI)と呼ばれるものです(“信頼”区間ではない;意味が異なるので注意)。先の真の値(alpha = 0.1, beta = 0.8)と比べると、いずれ真の値も95%CIに含まれていることが分かります。また、これらの推定結果と合わせて、psrfというものも表示されていますが、これはパラメーターの推定値が十分収束しているかどうかを判定するための値(Gelman-Rubin statistics、もしくはR-hatと呼ばれる)です。この値が1.1を上回る場合は推定値が十分に収束していないので、run.jagas関数のburnin、sample、thinの値を適宜調整し、信頼できる推定値を得る必要があります。

RからBUGSソフトウェアを動かす

BUGSファイルの記述で作った統計モデルを使い、JAGSによるパラメーター推定をしてみます。まずは、真のパラメーターが分かっている確率分布から、適当なテストデータセットを作ります。

###R code###
N.sample <- 100 # sample size
Area <- rpois(N.sample,100) # forest area
X <- scale(Area)  # forest area

alpha <- 0.1 # intercept
beta <- 0.8  # slope

y <- lambda <- NULL
for(i in 1:N.sample){
 lambda[i] <- exp(alpha + beta*X[i])
 y[i] <- rpois(1,lambda[i]) # Poisson error process
}

ここで作られた鳥の個体数のデータセットは、標準化された森林面積Xに対し、切片0.1(α)、傾き0.8(β)という関係性をもっています。このデータをプロットすると、100地点分の鳥個体数と森林面積の関係は右図のようになります。

このデータセットからパラメーターを逆推定し、真の値と比較することで、推定結果の妥当性を判定できます。データセットをJAGSへ渡し、計算を実行します。

###R code###
#①パッケージrunjagsの読み込み(JAGSを動かすための関数がある)
library(runjags)

#②JAGSに渡すためのデータセット
Djags <- list(N.sample = N.sample, y = y, X = X)

#③推定結果を知りたいパラメーター
para <- c("alpha","beta")

#④使う統計モデルの読み込み
m <- read.jagsfile("model_PoissonRegression.R")

#⑤数値計算のための初期値の設定
inits <- replicate(3,
                   list(.RNG.name="base::Wichmann-Hill",
                   .RNG.seed=round(runif(1,1,100))),
                   simplify = FALSE)

#⑥関数run.jagasでJAGSを動かし、計算結果をpostに格納
post <- run.jags(model = m$model, monitor = para, data = Djags,
                 n.chains = 3, inits=inits, method = "parallel",
                 burnin = 500, sample = 500, adapt = 100, thin = 3,
                 n.sims = 3)

①パッケージrunjagsの読み込み

install.packages(runjags)と打つことでインストールできます。このなかに、read.jagsfile, run.jagsという関数が入っています。

②JAGSに渡すためのデータセット

モデルを記述したBUGSファイルで使われているデータセットを入力します。ここで、モデルで使われているデータが入力されていないと、エラーがでて計算がストップします。なお、JAGSの場合(ほかの言語もほぼ同様)、データの型に以下のような制限があります。

  • データセットはlist形式(ただし、BUGSソフトウェアによって異なる)
  • 変数の中身はスカラー、ベクター、マトリクス、あるいはアレイ形式(X = 1, 2, 3, … nのような純粋な数字の羅列;変数名などが入るとエラーが出る)
  • 因子型の変数(“Yes”など)は認識されないため、すべて数字に変換しなければならない

③推定結果を知りたいパラメーター

推定結果を知りたいパラメーターを宣言します。ここに含まれていないパラメーターの推定結果は表示されません。なお、パラメーターはモデルファイル(model_PoissonRegression.R)で定義されているパラメーター名と一致する必要があります。

④使う統計モデルの読み込み

使うモデルファイル(model_PoissonRegression.R)を指定します。

⑤数値計算のための初期値の設定

JAGSのベイズ推定の過程では、MCMC(マルコフ連鎖モンテカルロ)と呼ばれる逐次的な推定手法が用いられています。パラメーターには初期値を与えると、指定した統計モデルのもとであるパラメーターのフィッティングを評価し(尤度)、次のステップではよりいいパラメーターの値を探す、という計算過程を繰り返し行います。ベイズの推定結果は、この(数値)計算過程の結果、パラメーターが収束する「値の範囲」を示したものになります。通常、このMCMCという計算は一度の推定で3つ行うのですが、それぞれ違う初期値からスタートする必要があります(違う初期値から同じ値に落ち着くかどうかを確かめたいため)。ここのコードは、その初期値が違う値となるよう設定する「おまじない」です。

⑥関数run.jagasでJAGSを動かし、計算結果をpostに格納

関数run.jagsが実際に「JAGSを動かせ」という指令を出しており、この関数が入力された時点でJAGSの計算がスタートします。それぞれの引数について説明します。

  • model: 実際に使う統計モデルを指定する。ここでは、「m」というところにモデルファイルが格納されているので、その中のm$modelを渡す。
  • monitor: 推定結果を確認したいパラメータを指定する。③で作ったparaを渡す。
  • data: 推定に使うデータセット。②で作ったDjagsを渡す
  • n.chains: MCMCの数を指定する。通常は3つで、変更する必要はほぼない。
  • inits: MCMCで使う初期値。⑤で作ったinitsを渡す。
  • method: 計算の際の詳細指定。いくつかあるが、parallelとすることでパソコンの複数のコアをフル活用し、計算が高速化される。
  • burnin: MCMCの最初の何試行分を捨てるかを指定する。ここで指定されたMCMCサンプルは推定結果に反映されない。全試行数の30~50%で指定することが多い。
  • sample: 最終的な推定値の計算に用いるMCMC1つ当たりのサンプル数。500とすると、500*3で1500サンプルが推定に用いられる。500~1000くらいが妥当。多いほど推定結果が収束しやすくなるが、その分計算に時間がかかる。
  • adapt: burninよりも先になされる試行計算。推定結果には反映されない。100~1000で設定しておけば問題ない。
  • thin: MCMCサンプルを「何個おき」に取ってくるかを指定する(その間のサンプルは「間引き」され、推定結果に反映されない)。MCMCの計算が自己相関(次のステップの値が前の値と類似する)するときには、間隔を広げたほうがいい。最低3は必要で、複雑なモデルになると1000以上の値を設定する必要がある場合もある。間隔をあけるほど全試行数は増えるので、計算時間は長くなる。
  • n.sims: 計算に使うCPUコア数の指定?のようなもの。chainの数に合わせる。

続きはこちら

BUGSファイルの記述

まずはBUGS言語による統計モデルの記述から説明を進めます。BUGSファイルは「統計モデルの構造」のみを記述するファイルです。ここでは、森林パッチiの鳥の個体数Yiを、森林面積Xiで説明するポアソン線形モデルを例に挙げて考えます。

Eq. 1)
Yi ~ Pois(λi)
log(λi) = α + β*Xi

上記のモデルは、RのGLMでは次のように記述されます。

###R code###
data <- read.csv(“data.csv”)
glm(y ~ X, data, family=poisson)

ここで気を付けなければならないのは、

  • GLMのようにRでパッケージ化されているものだと、eq.1で書かれているほとんど部分が自動化されており、リンク関数(log)などを自分で考える必要がない
  • 特に指定しなくとも、サンプルの繰り返しが“i”として認識されている

ことです。これら2点がBUGSとR言語では大きく異なります。このことを念頭におき、上記のモデルをBUGS言語で書いてみます。

###BUGS code###
model{
#Likelihood
 for(i in 1:N.sample){ # loop for i
   y[i] ~ dpois(lambda[i]) # poisson model
   log(lambda[i]) <- alpha + beta * X[i] # linear model structure
 }

#Prior
 alpha ~ dnorm(0, 0.001)
 beta ~ dnorm(0, 0.001)
}

model{}で囲われている範囲がモデル構造を指定するエリアになります({}の外にテキストを書くとエラーがでるので注意)。「~」は確率論的ノード(Stochastic node)と呼ばれ、左辺が確率変数である場合に使われます。確率論的ノードの右辺には、必ず確率分布の記述が必要です(上記ではdpoisが確率分布、Poisson distributionを指すもの)。「<-」は決定論的ノード(Deterministic node)と呼ばれ、=(イコール)と対応します。こうしてみると、BUGS言語ではEq.1をすべて自分で書き下す必要があると分かります。 Prior以降の記述は「事前分布」と呼ばれるものです。ベイズ統計の枠組みでは、すべてパラメーターに対し、事前分布を設定しなければなりません(事前分布を設定し忘れるとエラーが出る)。ここでは切片αと傾きβに対し、事前の情報は何も得られないので、「無情報事前分布」というものを設定しています(ほとんどケースはこれにあたる)。ここでは、平均0、分散1000の正規分布を事前分布として設定することで、切片αと傾きβは「かなり広い範囲の値をとる可能性がある」としています。「dnorm」が正規分布であることを示しており、一つ目の引数が「平均」、二つ目の引数が「精度(分散の逆数=1/1000)」を表します。ここまで書き終えたら、ファイルを.Rもしくは.txtファイルとして保存します(ここでは“model_PoissonRegression.R”として保存することにします)。続きはこちら

BUGSの使い方

最近ではRパッケージも実装され、容易にRからBUGSソフトウェアを使うことができます。その一連の流れは以下のようになっています。

  • 扱いたい統計モデルをBUGS言語で記述し、.Rもしくは.txtファイルとして保存
  • 解析に使うデータをRで準備(データ整形のしやすさはRのほうが優れているため)
  • RからBUGSソフトウェアにデータを渡し、指定した統計モデルのもとで計算を実行
  • BUGSからRに渡される計算結果を保存

この際、「どんな統計モデルを計算させるのか」という部分だけ、RではなくBUGS言語で記述する必要があります。したがって、統計モデルの実装の際には、RおよびBUGSスクリプトをそれぞれ別個に用意することになります。

ベイズ統計とBUGS

ベイズ統計、と聞くとなにやらすべてが新しいものに聞こえますが、扱うモデルの「構造」自体は大きく変わりません(線形回帰、ANOVAなど)。ただし、回帰係数のような“パラメーター”の扱いと推定手法が大きく異なっており、ベイズ統計の実装にはマルコフ連鎖モンテカルロのような逐次的な推定法が必要です。自前でアルゴリズムを構築することも可能ですが、そんな煩雑なプログラムを自分で書くというのは現実的な選択肢ではないでしょう。

そこで、BUGSソフトウェアの登場です。BUGSとはBayesian inference Using Gibbs Samplerの略で、ギブスサンプリングによるベイズ統計の実装を指します。このギブスサンプリングを実行してくれるのがWinBUGSをはじめとするBUGSソフトウェアで、Rから動かすことができます。そのため、昨今のベイズ統計を用いた研究では、RとBUGSソフトウェアの連携が一般的な手法になっています。なお、WinBUGSのほかにも、OpenBUGSとJAGS(Stanもベイズ実装のためのソフトですが、ギブスサンプリングではない)があります。WinBUGSは開発が終了していますが、OpenBUGSとJAGSは最新バージョンが随時ネット上でダウンロードできるようになっており、計算効率もそれらのほうがより良いものとなりつつあります。ここでは、BUGの中ではもっともクラッシュしにくいJAGS(Just Another Gibbs SamplerなのでBUGSというには語弊がありますが。。。)について解説します。

JAGSのダウンロードは次のページから可能です(JAGS)。ダウンロードはProgram files(デフォルト)もしくはCドライブ直下がお勧めです。