複数のパラメータ

Published on by anesth

Bayesian Biostatistics by E Lesaffre & AB Lawson, Chapter 4 学習ノート ???????????????????????
複数のパラメータをベイズ流ではどう扱うのかについての解説である。
たいへんハードな内容で、1日1章とはいかない。完全に理解するのに一ヶ月くらいかかりそう。
この際、慌てずにゆっくり行こう。一通りのプログラムをまずは動かして、何が課題なのかをざくっと理解することととする。

ALP.txtのデータ構造 315個のデータが収納されている。
artikel=0が250個、artikel=1が65個。alkfos[artikel==0]でartikel=0の250件を取り出しているが、残りのartikel=1が何か解らないが、ともかく250件のデータを分析する。

分散の大きさから、明らかに過分散データであるから、以下のように、正規分布に近似される100*√alp と整形している


頻度論アプローチで、N(μ, σ^2)に関して、観測値yの95%参照区間は、[μ-1.96σ, μ+1.96σ]
y_mean = 7.11で、s = 1.4なので、alpの尺度で95%参照区間は[104.45, 508.95]となるが、
この区間はμを推定するときのばらつきを無視している。yの標本分布を考慮すれば、y-y_mean〜N[0, σ^2(1+1/n)]であり、95%参照区間は[y_mean – 1.96σ√(1+1/n),y_mean + 1.96σ√(1+1/n)]となり、スケール化yは[4.43, 9.79]、もとの値では[104.33, 510.18]となる。

y=100/√alpとした同時分布p(μ,σ^2|y)を推測


と、美しい三次元立体図が描けた。
mugrid$muとmugrid$sigma2が、xy座標のパラメータで、z座標がヒストグラム頻度であろう。

mugridを排出するexpand.grid()関数とは何か？

expand.grid {base} R Documentation
Create a Data Frame from All Combinations of Factor Variables

Description

Create a data frame from all combinations of the supplied vectors or factors. See the description of the return value for precise details of the way this is done.

Usage

expand.grid(…, KEEP.OUT.ATTRS = TRUE, stringsAsFactors = TRUE)
Arguments

…
vectors, factors or a list containing these.
KEEP.OUT.ATTRS
a logical indicating the “out.attrs” attribute (see below) should be computed and returned.
stringsAsFactors
logical specifying if character vectors are converted to factors.
Value

A data frame containing one row for each combination of the supplied factors. The first factors vary fastest. The columns are labelled by the factors if these are supplied as named arguments or named components of a list. The row names are ‘automatic’.

Attribute “out.attrs” is a list which gives the dimension and dimnames for use by predict methods.

Note

Conversion to a factor is done with levels in the order they occur in the character vectors (and not alphabetically, as is most common when converting to factors).

References

Chambers, J. M. and Hastie, T. J. (1992) Statistical Models in S. Wadsworth & Brooks/Cole.

See Also

combn (package utils) for the generation of all combinations of n elements, taken m at a time.

Examples

require(utils)

expand.grid(height = seq(60, 80, 5), weight = seq(100, 300, 50),
sex = c(“Male”,”Female”))

x <- seq(0, 10, length.out = 100) y <- seq(-1, 1, length.out = 20) d1 <- expand.grid(x = x, y = y) d2 <- expand.grid(x = x, y = y, KEEP.OUT.ATTRS = FALSE) object.size(d1) - object.size(d2) ##-> 5992 or 8832 (on 32- / 64-bit platform)
ーーーーーーーーーーーーーーーーーーーーーー

は何をしているのか？

サンプルデータの作成：関数 expand.grid()

関数 expand.grid(引数1=ベクトル1，引数2=ベクトル2，・・・) を用いることで，引数に指定した因子全ての組み合わせから成るデータフレームを作ることができる．

expand.grid関数を使って、引数のベクトルの全組み合わせを作ることができます。例えば、こんな風に作れます。
Rコンソール
> expand.grid(1:2,1:3)
Var1 Var2
1 1 1
2 2 1
3 1 2
4 2 2
5 1 3
6 2 3

expand.grid()の中身を見てみる。
mu=do.breaks(c(unl.mu,upl.mu),nmu),
sigma2=do.breaks(c(unl.sigma2,upl.sigma2),nsigma2)

平均が7.1なので、7.1 ± 0.3をunl.mu, up1.muとしているようだ。
分散は1.86、標準偏差は分散のルートで1.36。
unl.sigma2 <- 1.2, upl.sigma2 <- 2.5がどこから来たのかわからない。

次のパッケージでは、

mu <- seq(6.2,7.5,length=100) mus <- (mu-7.1)/(σ/sqrt(n)) pmu <- dt(mus,n-1)*sqrt(n)/σ sigma2 <- seq(1.2,4,length=100) sigma2s <- (n-1)*var/σ2 psigma2 <- dchisq(sigma2s,n-1)*(n-1)*var/(σ2^2)

μの事後分布の位置パラメータはμ=7.11で、分散σ_mean^2 = 0.0075
95%信用区間は[6.94, 7.28]
σ^2の事後分布の要約値はσ_mean^2=1.88、最頻値1.85、中央値1.87で、σ_mena_σ=0.029、σ^2に値する95%信用区間は[1.58, 2.24]で、95%HPD区間は[1.56,2.22]

将来の観測値の分布はt249(7.11, 1.37)分布となる。したがって、yの95％正常範囲は[4.41,9.80]であり、alpの95%正常範囲は[104.1, 513.2]となる。

次に、後ろ向きの事前研究のデータalpのデータのalkfos[artikel==1]の65個のデータを使う。
共役事前分布はN-Inv-χ^2(5.25, 65, 64, 2.76)
事後分布はN-Inv-χ^2(6.72, 315, 314, 2.61)
事後平均6.72は、事前平均5.25と標本平均7.11の重み付き平均である。μの事後分散は0.0083であり、以前の無情報事前分布から求めた事後分散0.0075より大きい。
μの95%信用区間は[6.54, 6.90]
σ^2の事後平均は2.62で、事後分散は0.044、95%信用区間は[2.24, 3.07]で、95%HPD区間[2.23, 3.04]
将来の観測値の分布はt314(6.72, 1.62)分布であり、yの95%信用区間は[3.54, 9.91]で、対応するalp 95%正常範囲は[101.9, 796.7]である。

4.4 多変量分布
青年調査：喫煙と飲酒の関連


Kaldor’s et al. のケース・コントロール研究：例?.9の再解析

4.5 ベイズ推測の頻度論的な性質

4.6 事後分布からのサンプリング：合成法

4.7 ベイズ流線形回帰モデル

4.8 ベイズ流一般化線形モデル