2021.03.19
丹後俊郎『統計モデル入門』(朝倉書店)は10年前に癌の治療成績の医療統計を始めた時に買った本であるが、当時はまだ僕には難しかった。同じ著者の『医学への統計学』をもっぱら頼りにして、R で計算していた。現在、新型コロナ感染症の解析を試みる中で、いろいろと新しい統計手法が使われているので、それらを理解する為につまみ食いして読んでいる。(自分でプログラムまで組むつもりはないので、、。)実用に即して本質的な説明がしてあって判りやすい。随時まとめを追記していく予定である。おや、と思ったのは第2章で、食中毒の発症日の分布のデータだけから、感染日まで特定できるという。確かに食中毒の場合には感染が一回しかないのだから、潜伏期間が対数正規分布と仮定さえすれば、発症日の分布だけで、食中毒の感染日が推定できるのである。。。

● 第3章のブートストラップは最近良く出てくる概念である。読んでいて、Bias というのが良く判らなかった。これはブートストラップサンプルでのパラメータ期待値と実際のサンプル(これをそのまま分布 F1 と見做す)でのパラメータ推定値との差の事である。ブートストラップというのは、観測データからランダムに重複を許して再サンプルすることで、これをそのまま本当の母集団に見立てるものだろうと、僕は勝手に思っていたから、単に前者の推定値を採用すればよいだけで、どうして差が問題になるのだろう、と思った次第である。これは僕の勘違いであった。

・・・現実に存在するデータ X は何らかの分布 F0 に従う確率変数であると考える。その確率分布 F0 を規定するパラメータ θ を推定する。このとき、データ X 自身は分布 F1 を構成しているから、その分布におけるパラメータ θ(これをθ~ と書く)が得られる。これが F0 のパラメータ θ とどれくらいずれているか、がバイアスである。F0 は知られていないから、これは直接計算することができない。そこで、そのバイアスを推定するために、今度は F1 を元の分布と見做して、そこから F2 を抽出することで模擬的に状況を作り出す。それは X 自身からランダムに選択して多数の疑似分布を作り出し、それらから得られるパラメータ θ の期待値がどれ位 θ~ からずれているか、ということで推定するのである。(これは現代の計算機ではいくらでも乱数を発生して容易に行える。)つまり、

    バイアス≒E(θ(F2)) - θ~ :E() は期待値を意味する。

である。そうすると、このバイアスの分だけ、θ~ が本当の θ(F0)からずれていると推定されるから、

    θ~ ≒ θ(F0) + E(θ(F2)) - θ~
となり、
    θ(F0) ≒ 2θ~ - E(θ(F2)) が得られる。

・・・ブートストラップ分布でずれたのとは逆方向に補正することになる。実例を見るとへえーッと驚くかもしれない。
11人の患者のGPT(血液検査項目)のデータがあって、その中央値は θ~=129である。具体的には、
{16, 41, 57, 76, 124, 129, 193, 215, 280, 363, 914}
ここからランダムに選択して11個のデータを作り出して、その中央値を求める、ということを100回、200回、、、と繰り返して、その中央値の期待値を求める。これが、E(θ(F2))で、149.8に収束する(実は解析的にも計算できる)。バイアス推定値は 20.8 ということになる。これが F0 と F1 の間のバイアスの推定値と見做されるから、F0 における中央値は 108.2 と推定されるのである。サンプルを観察すると、F0 の分布は小さい値の頻度が高くて、大きな値の頻度が小さい、つまり大きな値の方に裾を引いていると思われる。そういう分布から少数のサンプルを採択した場合には、その少数サンプルの中央値が本当の中央値よりも高めになる、ということを意味していて、その傾向(バイアス)がブートストラップによってシミュレーションされている、ということだろうと思う。
・なお、平均値や(不偏)分散についてはバイアスは 0 である。

・・・信頼区間についても同じ論理である。θ が F1 から計算された θ~ からどれくらい外れている可能性(確率)があるかを知りたい。
    P(θ~ - t0 < θ(F0) < θ~ + t0) = 0.95。
しかし、本当の分布 F0 も本当のパラメータ θ も判らない。そこで、F0 の代用として F1 を採用し、F1 の代用として F2 を採用することで、これを模擬的に再現して、t0 を推定する。
    P(θ(F2) - t0 < θ~ < θ(F2) + t0) = 0.95。
θ(F2) は乱数を発生させて、X からランダムに選択したデータを作り出せばいくらでも作り出せるから、その分布関数が判る。式を変形して、
    P(θ~ - t0 < θ(F2) < θ~ + t0) = 0.95。
つまり、θ~ の上下に t0 だけ区間を拡げて、その中に、乱数で発生させた θ(F2) の内、95% が入るように、t0 を選ぶ。これをパーセンタイル C.I 法というらしい。

・・・しかし、Bias がある場合はどうなるのだろう?
θ(F2) を 1000回繰り返し、小さい方から順に並べたとする。θ(1)~θ(1000) である。
    P(θ(25) - θ~ < θ(F2) - θ~ < θ(975) - θ~)=0.95
で θ(F2) - θ~ を θ~ - θ(F0) と見做して、
    2θ~ - θ(975) < θ < 2θ~ - θ(25)
となる。これは確かに Bias 補正と整合する。これは標準 C.I と呼ばれているらしい。これらの C.I.の算出方法を更に改良した方法がいくつかある。。。

・・・Bootstrap サンプルというのは、母集団をシミュレートするものだろうと思っていたのだが、どうもそれではないようである。むしろ可能性のある標本集団をシミュレートするものである。現実に観測されている一つの標本を仮に母集団と見做したとき、そこから標本を観測したとしたら、どういう標本が得られるか、ということである。それを計算機上で仮想的に作り出すことで、標本を得る時の Bias やバラツキを推定し、その Bias やバラツキが本当の母集団と観測された標本との間にあると推定するのである。何だか騙されたような気分になるが、漸近的には正しいことが証明されているようである。

● 第9章がベイズ統計である。最初に頻度統計の概念が説明してある。
データから一つ存在すると想定される母集団を規定するパラメータを求める、という立場である。
例。。。
患者 n 人(i=1,n) と測定 r 回(j=1,r) のデータ(一元配置) Xij を

    Xij = μ + αi + εij  :εij は誤差。

というモデルで考える時、μ と αi が母数である。これを決める時、固定効果モデルというが、

    αi ~ N(0,τ)

として、αi そのものではなく、分散 τ を決めようとすれば、これは変量効果モデルである。
(新型コロナでの無症状者数分布のメタアナリシスで使われていた。)
更に、測定の時間(tj)依存性と個体差を考慮して、

    Xij = (μα + αi) + (μβ + βi)tj + εij

とすることもできる。要するに、何に興味があるか、によってモデルが選択される。

・・・これに対して、ベイズ統計では、母集団を規定するパラメータ θ そのものが確率変数となる。これは何も分布を規定するパラメータでなくてもよい。欠損値の推定とか測定誤差を除いた数値とか。その確率分布に事前確率 p(θ) を仮定する。勿論観測値 X も確率変数であり、それは、θ という条件の元での確率分布 f(x|θ) に従う。この形は知られているとする。そうすると、θ であり、且つ x である、という同時確率が定義できる。これは x であり、且つ θ である、とも読める。つまり、

    p(θ,x) = p(θ)f(x|θ) = f(x)p(θ|x)

ここで、f(x) は x となる全確率、つまり、

    f(x) = ∫ p(θ)f(x|θ)dθ

である。従って、事後分布は、

    p(θ|x) = p(θ,x)/f(x)

となる。分母は θ を含んでいないから、分布の相対値だけを問題にするときには単なる定数因子である。

・・・実際上問題になるのは、事前分布の選択方法である。社会学では主観的であっても構わないかもしれないが、医学ではそういう訳にも行かない。客観性が要求される。無情報事前分布とempirical Bayes 法がある。

・・・無情報の場合。パラメータ範囲が正負に亘る場合、正規分布(以下 N(μ,τ) と表記する)が便利である。正に限られる場合には逆ガンマ分布が便利である。

  p(1/θ)=Ga(a,a)、a は0.001程度に採れば充分「客観的」である。なお、

  f(x|a,b)=b^a/Γ(a) x^(a-1) exp(-bx)。

これらの分布関数の便利な処は、事前分布と事後分布が同じ分布属になることである。

・・・事後分布(パラメータの分布)が得られたとして、最適パラメータをどうやって決めるか?

・・・損失関数 L(t,θ) を考える。その p(θ|x) での期待値を最小化するような x を最適な θ とする。

・L=(t-θ)^2 とすると、この場合 t はパラメータの期待値となる。これが標準的な考え方らしい。
・L=|t-θ| とすると、t はメディアンとなる、
・L=0 for t=θ、1 for t≠θ とすれば、モード(最頻値)で、これは最尤法と同じである。
つまり、ベイズ法は最尤法をその一部として含む。こういう処に違いがあった。。。

・・・階層構造を持つベイズモデルが応用上重要である。これはグラフで表現できる。未知パラメータを○で囲む。データは□で囲む。依存関係が→。実線は確率的依存性で、点線は関数的依存性。
(例)
平均値 μ と分散 1/τ を持つ未知分布からの確率変数のサンプル x から平均値と分散をベイズ推定する問題。
ただ、この例では正規分布を想定している。
 τ固定で、μ の事前分布は分散 1/τ' の正規分布とする。xi が同時に実現しているから確率の積となって、

  p(μ,τ,x)∝exp(-τ'μ^2/2)Πexp(-τ(xi-μ)^2/2)
          (Π は exp の中に入ると ∑になる、τ は定数、以下導出略)
         =N(n<x>τ/(τ'+nτ),1/(τ'+nτ))→N(<x>,1/nτ) for n→∞:<x>は観測データの平均値である。

  μ固定で、

  p(μ,τ,x)∝τ^(a-1) exp(-aτ)Πτ^(1/2) exp(-τ(xi-μ)^2/2)
          (この場合τが変数なので、規格化条件が変化する、以下導出略)
            ∝Ga(a+n/2, a+Σ(xi-μ)^2/2)

斯くして、事後分布も同じ分布属になるのである。

  μの期待値=n<x>τ/(τ'+nτ)
 τの期待値=(a+n/2)/{a+Σ(xi-μ)^2/2}
 これらを連立させて解くことで、ベイズ推定値が得られる。

● 第10章:マルコフ連鎖モンテカルロ(MCMC)法
・・・ベイズ推定には、上記のように簡単な場合を除いて、多重積分が必要となるので、解析的には難しくなる。
例。
・・・毒物の量(の対数) xk のそれぞれに対して、実験動物 nk 匹を用意して、その死亡数 mk を観測する。このデータに対して、毒物量に依存した死亡率 θk を近似するロジスティック回帰式を求める。死亡数は二項分布の確率変数と考える。

    mk ~ Binom(θk,nk)
    log(θk/(1-θk)) = α + βxk + εk

として、パラメータ α、β を求めればよい。
事前分布として p(α,β) を仮定する。βについての事後分布は

    p(β|x,m,n) = ∫p(α,β)f(m|x,n,α,β)dα/∫∫p(α,β)f(m|x,n,α,β)dαdβ

ここで、f は その死亡率パラメータ θ が α、β、xk と回帰式から推定される二項分布である。
必要なのは、β の期待値(ベイズ推定値)である。

    <β> = ∫βp(β|x,m,n)dβ

・・・一般的に、確率分布関数 π(x) による変数 f(x) の期待値 Eπ(f(X))=∫f(x)π(x)dx を計算するのに、
π(x) を模擬するようなランダムなサンプル xi を使って、

    Eπ(f(X))≒(1/N)∑f(xi)

と出来れば良い(Monte-Carlo法)。
ある分布に収束するような乱数の系列を効率的に(早い収束)作り出す方法として、マルコフ連鎖がある。

    x(i+1) ~ p(x|x(i))

つまり、一つ前の乱数にのみ依存する(マルコフ連鎖)確率分布に従う。
よく例に出されるのが、この確率分布として N(0.5x(i),1.0) とする場合である。つまり前回の乱数の半分を平均値とする正規分布に従うとする。これは平均が 0、分散が 4/3 の正規分布に収束する。

・・・Metropolis-Hastings の方法 MH法
  任意のサンプラー q(・|・) を設定して、まず

    y ~ q(x|x(i))

とする。

    α(x(i),y) = min(1, π(y)q(x(i)|y)/π(x(i))q(y|x(i)))

とする。つまり、x(i)→y という継起事象と y→x(i) という逆の継起事象の生成確率(積確率)を比較している。
順の方が大きい時は α=1 として、小さい確率の方に移動する時はその比を α として採用する。
そこで、新たな x(i+1) = y として更新する確率を α とする。つまり連続事象が大きければ全て更新する。
連続事象の確率が小さい場合には、α だけの確率で y へと更新する。
  (統計力学のシミュレーションでよく使われるのは、q が酔歩運動の場合である。π が熱平衡分布である。)

・・・いろいろなサンプラー
・対照サンプラー q(y|x) = q(x|y) とりわけ、N(X,σ^2)、random walk q(y|x) = q(|y-x|)。
・独立サンプラー q(y|x) = q(y)

・・・Single-component MH 法
・多変数(=パラメータ)の場合には、順に行う方法が使われる  Single-component MH 法
一つ一つの変数 xj の遷移に際しては、他の変数は現在の状態に固定するから、その固定した条件での確率分布が必要となる。
これをフル条件付き分布と呼ぶ。π(xj(i)|x-j(i)) : -j は j以外の変数を表す。
つまり他の変数で積分する必要がある。

・・・Gibbs sampling:現在は、殆どこれが使われているらしい。
・この π(xj(i)|x-j(i)) をそのまま qj(yj|xj(i),x-j(i)) として採用したものが Gibbs sampling と呼ばれている。
これは結果的には独立サンプラーであり、且つ α=1 となるので、常に採択される。
実際上複雑な確率分布を持つ乱数を発生させるのは容易でないから、rejection sampling が使われる。
π よりも 大きな関数 Π に比例した乱数を発生させて、独立に一様分布乱数によって、その中から π/Π の確率で採用する。
Π の選択方法として、代表点から secant method によって折れ線分布 Π を作る方法が adaptive rejection method である。

  <目次へ>       <一つ前へ>     <次へ>