生態学分野においてモデル平均は評判が悪い気がするが、(例：http://hosho.ees.hokudai.ac.jp/~kubo/ce/FaqModelSelection.html#toc9)
「データに条件づけられた偏回帰係数の分布を、各モデルの事後分布とモデルの事後確率で重みづける」ということ自体はベイズの定理だけ導出できる([1,2])。しかし、Rの"MuMIn"パッケージで計算される「重み付き推定値に基づくp-value」には、重大な数学的欠陥がある。以下で、その問題点を指摘したい。
目的変数をy(y1,y2,y3,...yn)、説明変数をX、偏回帰係数βとして複数の回帰モデル(M1,M2,M3,...Mi)を考えてみよう。
モデル平均を用いたβの事後分布は、
p(β|Data)=Σp(βi|Data)*p(Mi|Data)
ここで、p(βi|Data)はモデルiにおけるβの事後分布、p(Mi|Data)は、モデルiの事後確率である。
=================参考===============================
ただし、ふつうモデルの事後確率は計算上の困難を伴うため、AICやBICを用いた近似が有用である（*注1）。
具体的には、
p(Mi|Data)=p(Data|Mi)*p(Mi)/ Σp(Data|Mi)*p(Mi)
において、
p(Data|Mi)をexp(ΔAICi/2)かexp(ΔBICi/2)で近似する。
ここで、ΔAICiは候補モデル群の中で最も小さなAICを記録したモデルのAICと、モデルiのAICの差分である。（BICを用いた場合についても同様）
この方法でモデル事後確率を近似したのが、いわゆるAkaike weightやSchwartz weightである。
=====================================================
さてモデルの事後確率で重み付けされたβの事後分布は、事後分布が正規分布であると仮定すると（すなわち、モデルが正則であると仮定すると）以下の形でかける
p(β|Data)=
p(M1|Data)*Normal(μ1,σ1)+p(M2|Data)*Normal(μ2,σ2)+p(M3|Data)*Normal(μ3,σ3)+・・・+p(Mi|Data)*Normal(μi,σi)
p(Mi|Data)がモデルiの事後確率で、Normal(μi,σi)が、モデルiのもとで推定された事後分布である。
これは、いわゆる混合正規分布というものである。混合正規分布は、混合比パラメータと、各正規分布の平均・分散に応じて非常に劇的な形状の変化を見せる。
例えば以下は、混合させる正規分布を二つだけに絞り、且つsd＝1で一定と他場合に限定して、混合比と平均パラメータを適当に与えたものである。
単峰性の普通の正規分布で混合正規分布を近似するには、かなり危険が伴うように思われる。
ところが、Rパッケージの"MuMIn"で出力されるp-valueは、この混合正規分布を重み付き推定量と、その標準誤差を使って正規分布で近似できる、という前提に基づいて計算されている。Rの計算例を見ながら、このことを確認してみよう。
いま、ある重回帰モデルを"a"としてobjectに格納する。そのサブモデルについて、事後モデル確率の近似計算と、偏回帰係数の重み付き平均を計算する。
#例：
a<-lm(y~x1+x2+x3+x4+x5,data=data)
#これは線形回帰モデルでなくとも、(MuMInによるモデル平均に対応しているようなクラスのモデルであれば)どんな種類のモデルでも良い。例えば、glm()、glm.nb()、zeroinfl()、glmer()、glmmML()、MCMCglmm()・・・
#パッケージMuMInの呼び出し
library(MuMIn)
#以下のオプションを設定しないと、計算が回らない
options(na.action = "na.fail")
#モデルaのサブモデルすべてについてAIC近似に基づくモデル平均を行い、object”avg.model”に格納する。
avg.model<-model.avg(get.models(dredge(a,rank="AIC"), subset = TRUE))
#結果の呼び出し
summary(avg.model)
各モデルがAIC順に羅列された下に、各偏回帰係数の重み付き推定値と、SE、調整SE,z-valueそしてp-valueが記載される。
任意の回帰係数について、
推定値/調整SE
を計算すると、右側のz-valueと同じになることが確認できるはずだ。
さらにRに標準装備された正規分布の累積密度関数を用いて
(1-pnorm(推定値/調整SE))*2
とすると、表の再右に記載されたPr(>|z|)に一致することが確認できる。
*注1
モデルの事後確率を計算する上で、本来はモデルの事前確率と偏回帰係数βの事前分布が必要である。
ここでは、すべてのモデルの事前確率が等しいと仮定した。
βの事前分布についてであるが、AICを用いた場合は、パラメータの尤度関数と同じ関数を事前分布に採用したとみなすことができる。
BICの場合は、「単位情報事前分布」を採用したとみなすことができる。
引用：
[1]Hoeting, J. A., Madigan, D., Raftery, A. E., & Volinsky, C. T. (1999). Bayesian model averaging: a tutorial. Statistical science, 382-401.
[2] Burnham, K. P., & Anderson, D. R. (2004). Multimodel inference: understanding AIC and BIC in model selection. Sociological methods & research, 33(2), 261-304.