ニューファンドランドのグレートアイランドのツノメドリpuffinの繁殖の成功率を調べたデータ。
38羽の鳥について、営巣頻度NEST,草類被覆GRASS, 平均土壌深度SOIL, 崖傾斜ANGLE, 崖縁からの距離DISTANCEを観測。
NESTi = β0 + β1(DISTANCEi – DISTANCE平均）+ ei
eiは正規分布N(0, σ)からのランダム標本。回帰係数ベクトルβの事前推定を（8,0)とする。8は、平均営巣頻度からきている。
> data(puffin)
>
> puffin
Nest Grass Soil Angle Distance
1 16 45 39.2 38 3
2 15 65 47.0 36 12
3 10 40 24.3 14 18
4 7 20 30.0 16 21
5 11 40 47.6 6 27
6 7 80 47.6 9 36
7 4 80 45.6 7 39
8 0 15 27.8 8 45
9 0 0 41.9 8 54
10 0 20 36.8 5 60
11 15 40 34.9 31 3
12 21 60 45.2 37 12
13 12 95 32.9 24 18
14 8 50 26.6 11 24
15 9 80 32.7 10 30
16 6 80 38.1 5 36
17 16 60 37.1 35 6
18 25 60 47.1 35 12
19 13 85 34.0 23 18
20 13 90 43.6 12 21
21 11 20 30.8 9 27
22 3 85 34.6 6 33
23 0 30 37.7 8 42
24 0 75 45.5 5 48
25 0 15 51.4 8 54
26 18 40 32.1 36 6
27 19 40 35.4 37 9
28 8 90 30.2 11 18
29 12 80 33.9 9 24
30 10 80 40.2 11 30
31 3 75 33.5 7 36
32 0 65 40.3 10 42
33 0 60 31.4 5 39
34 0 70 32.7 2 48
35 0 35 38.1 8 51
36 0 80 40.3 12 45
37 0 50 43.1 13 51
38 0 50 42.0 3 57
> X=cbind(1, puffin$Distance – mean(puffin$Distance))
> c.prior=c(0.1,0.5,5,2)
> fit=vector(“list”,4)
> for (j in 1:4)
+ {
+ prior=list(b0=c(8,0), c0=c.prior[j])
+ fit[[j]]=blinreg(puffin$Nest, X, 1000, prior)
+ }
> BETA=NULL
> for (j in 1:4)
+ {
+ s=data.frame(Prior=paste(“c =”,as.character(c.prior[j])),
+ beta0=fit[[j]]$beta[,1],beta1=fit[[j]]$beta[,2])
+ BETA=rbind(BETA,s)
+ }
> library(lattice)
> with(BETA,xyplot(beta1~beta0|Prior,type=c(“p”,”g”),col=”black”))
ここで４通りの値を与えたｃは平均β1と、分散共分散行列V1の計算式に含まれる係数であるが、事前分布のデータの情報量を反映させる定数で、事前の確信に重きをおくなら小さな値、逆に大きな値を選ぶと、（β、σ2^2)に標準的な無情報事前分布を選んだ場合と同じ効果となる。
blinreg()関数は、priorオプh村を指定すると、Zellnerのg事前密度を使って回帰モデルからシミュレーション標本を出力する。
Zellnerのg事前密度クラスの特徴は、回帰分析の問題で最適なモデルを選択するために使用できる。
反応変数yにk個の予測変数の候補があるとする。すべての予測変数を含むフルモデルをβと表す。
このβに事前の推測β0=0によるg事前分布と大きめのc (c=100)を割り当てる。すなわちβについては、曖昧な事前情報しか無いということである。その上で、βPが予測変数の部分集合Pを含む回帰モデルだとすると、βPに事前推測が0で同じcである同一のg事前分布を割り当てる。
異なる回帰モデル間の比較は、事前予測密度のの計算で行う。反応変数のサンプリング密度にf(y|β, σ^2)を割り当て、パラメータのベクトル(β, σ^2)に事前密度g(β,σ^2)を割り当てると、yの事前予測密度は以下の積分となる。
m(y)=∫f(y|β, σ^2)g(β, σ^2)dβdσ^2
ここで共変量GRASSとSOILを含む特別なモデルの事前予測密度を計算したいとする。
はじめに反応ベクトルyとデザイン行列Xを含むリストdataを定義する。
> data=list(y=puffin$Nest, X=cbind(1,puffin$Grass,puffin$Soil))
βの適当な初期値はlaplace()関数で最小二乗法による回帰あてはめによって求める。
reg.gprior.post()関数はg事後密度による回帰モデルの対数事後密度を計算する。
ここで、対数事後密度は、対数尤度と対数事前密度の合計。
> prior=list(b0=c(0,0,0), c0=100)
> beta.start=with(puffin,lm(Nest~Grass+Soil)$coef)
> laplace(reg.gprior.post,c(beta.start,0),list(data=data,prior=prior))$int
[1] -136.3957
bayes.model.selection()関数は、上のアルゴリズムを使って、すべてのk個の予測変数の組み合わせのモデルの予測密度を計算する。
> X=puffin[,-1]; y=puffin$Nest; c=100
> bayes.model.selection(y,X,c,constant=FALSE)
$mod.prob
Grass Soil Angle Distance log.m Prob
1 FALSE FALSE FALSE FALSE -132.18 0.00000
2 TRUE FALSE FALSE FALSE -134.05 0.00000
3 FALSE TRUE FALSE FALSE -134.51 0.00000
4 TRUE TRUE FALSE FALSE -136.40 0.00000
5 FALSE FALSE TRUE FALSE -112.67 0.00000
6 TRUE FALSE TRUE FALSE -113.18 0.00000
7 FALSE TRUE TRUE FALSE -114.96 0.00000
8 TRUE TRUE TRUE FALSE -115.40 0.00000
9 FALSE FALSE FALSE TRUE -103.30 0.03500
10 TRUE FALSE FALSE TRUE -105.57 0.00360
11 FALSE TRUE FALSE TRUE -100.37 0.65065
12 TRUE TRUE FALSE TRUE -102.35 0.08992
13 FALSE FALSE TRUE TRUE -102.81 0.05682
14 TRUE FALSE TRUE TRUE -105.09 0.00581
15 FALSE TRUE TRUE TRUE -101.88 0.14386
16 TRUE TRUE TRUE TRUE -104.19 0.01434
$converge
[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[16] TRUE
最も適切なモデルは{SOIL, DISTANCE}で対数周辺密度と事後確率はそれぞれお-100.37と0.65065となる。