2021.03.22
（駄文）
丹後俊郎『メタ・アナリシス入門』（朝倉書店）を借りてきていて、序文を読んで、ああそうか、と思ったのだが、統計モデルというのは、物理モデルとは随分違う。物理モデルの場合には本来的にそれよりも下位にある基本法則から説明できる、つまり何らかの現実的近似の元でそのモデルが成り立つ、という「メカニズム」が想定されているのだが、統計モデルにはそういう暗黙の制限が無い。一元論から切り離されている。プラグマティズムとも言える。「黒い猫だろうと白い猫だろうと鼠を捕る猫なら良い」のである。哲学だってそういうものである。生きる為の技術、あるいは方便に過ぎない。

・・・こういう統計モデル的発想は、企業活動においては容認せざるを得ないものであるし、僕が企業で学んだことも結局はそういうことであった。そこで原理原則に拘っていては相手にされない。物理法則や化学法則から見てなかなか信じがたいようなモデルであっても、現実のデータがその傾向を示しているならば、当面それを信じて、次の対策を講じるべきである。そうすれば新しいデータが得られる。その時点でモデルを更新すればよいのである。疫病や公害などでも同様である。そこで「病理的メカニズムが判らないから認められない」という意見が通ってしまえば、対策が遅れて水俣病のような悲惨なことになってしまう。つまり、「メカニズム」は判らない場合が殆どなのだから、データが統計的に傾向を示しているならば、それを否定してはならないのである。これが疫学の原理である。だからこそ、その傾向をどうやって引き出すか、とか、その傾向がどの程度真実性があるか、ということを「メカニズムが想定できない状況において」「データだけから」判断しなくてはならない。それが統計学の目的である。

・・・物理モデルの究極は超ひも理論 と目されているが、まだ実証されていない。今の処場の量子論に基づく標準理論までは確実である。統計モデルの方は究極が見えない。それは計算機の能力次第でどこまでも拡がる。むしろそちらの方で人間の知性が試されているのかもしれない。哲学者達は右往左往するばかりである。

（本文）
新型コロナウィルスの無症状者陽性者数比率とか、無症状者陽性者中での不顕性感染者数比率とかが必要になって、いろいろな文献を探し、またそれらの文献を統合して解析した文献を探した。この統合のことをメタ・アナリシスと呼ぶらしいのだが、そのやり方として殆どの場合に、変量モデル（random effect model）が使われている。具体的には DerSimonian and Laird モデルが使われているようなので、その計算方法をネットで調べて、エクセル上で計算をチェックして、やっと結果を利用することができたのだが、そのやり方がどうやって導かれたのがが判らないままである。困ったことに原論文（1986年）を読んでも判らない。これはやはり適切な教科書を読むしかない、ということで、丹後俊郎『メタ・アナリシス入門』を借りてきたのである。中身の殆ど（１～８章）は具体例と具体的にどの方法をどう使えばよいのか、の解説である。最後の９章で、統計学的にそれらの方法が導かれた論理を説明している。

・・・まずは母数モデル（fixed effect model） である。ある医学的手法の効果をさまざまな施設で試みた結果の報告を統合する。それらの結果を効果として θi で表すが、実際の報告はその（本当の隠された効果、一般的には確率分布のパラメータの）期待値（平均値）<θi> である。i は報告数 K までの自然数である。この効果の表現方法としては治験数に対する効果ありの数の比率でも良いし、差でも良いし、効果あり対効果無しの比率でも良いし、その対数でも良い。他にもいろいろな表現方法があるが、それらを変換して漸近的に正規分布となるような変数 f(θ) にしておくとする。以下式に全て f() が入ってくるので、ここでは省略して、単に θ とする。（なお、僕の場合、θ は無症状陽性者中の不顕性感染者数比率であったから、二項分布を正規分布と近似する必要があった。）つまり、推定値 <θi> にはその漸近分散、尤度分布の分散、 si^2 がある。

（注）二項分布の場合----------------------
患者 n 人中 k 人に効果があったとすると、<θ>＝k/n である。これは推定された二項分布パラメータ p である。二項分布において効果あり人数 k の分散は p(1-p)n となるから、二項分布を正規分布として近似すると、k ～ N(pn, p(1-p)n) であり、具体的には ～ exp(-(k-pn)^2/2p(1-p)n) である。今、観測されているのは k であり、それを前提として p を推定するときの尤度分布と見做すから、指数内分母の分散部分は観測値から p≒k/n と近似して、

    p ～ exp(-(p-k/n)^2/(2k(n-k)/n^3))
となる。従って、

    s^2＝k(n-k)/n^3
とすればよい。二項分布そのものの推定分散 k(n-k)/n ではない事に注意すべきである。 
このような変換が必要なので、f(θ) と表記してある。
------------------（注終り）

ここで便宜上（式が簡単になるように）

    wi = 1/si^2
と定義しておく。
帰無仮説は、効果が無い、つまり母集団においてそのパラメータとしての効果期待値 θi = θ = 0 である。この確率が小さければ、治療効果がある、という結論になる。これはそれぞれの試験施設において偏りが無いということでもある。
対数尤度は（正規分布だから）

     L(θ) ∝ Q = ∑i wi(<θi> - θ)^2
となる。< > はサンプル平均値。
これを最大化し（θで微分して 0)、その先鋭度を計算（θ で２回微分して半分）すれば、推定値（~で表す）と分散と信頼区間が得られる。結局の処、分散の平方根の逆数（∝サンプルサイズ）が重み因子になっている。

     θ~ = ∑i <θi>wi/(∑i wi)
    分散 ＝ 1/∑i wi
     95% c.i → ±1.96/√∑i wi
となる。結果的には全ての測定値を同等と見做して、その平均と分散を計算していることになる。ところで、

     Q = ∑i wi(<θi> - θ~)^2 + ∑i wi(θ~ - θ)^2
と分解できて、第１項 Q1 は自由度 K-1 の χ2 分布に従い、第２項 Q2 は自由度 1 の χ2 分布に従う。
つまり、それぞれの量が χ2 分布のどこにあるかを見れば、第１項によって、試験施設 K 個の間の均質性が検定されて、
第２項によって、統合された効果の有無（これが一番の目的）が検定される。
均質性が疑われる（Ｑ1 が χ2 分布の端の方にある）と、そもそもこの母数モデル自身が成立しなくなるのでモデルを変えなくてはならない。

・・・そこで、変量モデル（Random Effect Model）が必要になる。
ここでは、試験施設による偏りがあり、その為に本当の効果 θi がばらついているとする。その分散 τ^2 （τ はギリシャ文字 タウ）を考えると、正規分布の積となるから、

     <θi> ～ N(θi,si^2 + τ^2)
となる。θi が θ の周辺で分散 τ^2 の正規分布をして、その θi の周辺で <θi> が si^2 の正規分布をする、と考えると、

     <θi> ～ exp{-(<θi> - θi)^2/2si^2} exp{-(θi - θ)^2/2τ^2}
          = exp[{-<θi>^2/si^2 - θ^2/τ^2 - (1/τ^2+1/si^2)θi^2 + 2(<θi>/si^2+θ/τ^2)θi}/2]
          = exp[{-<θi>^2/si^2 - θ^2/τ^2 
                     - (1/τ^2+1/si^2){θi - (<θi>/si^2+θ/τ^2)/(1/τ^2+1/si^2)}^2
                     + (<θi>/si^2+θ/τ^2)^2/(1/τ^2+1/si^2)}]/2

となり、これを θi で積分すれば、単なる確率変数和の正規分布（平均も分散も和になる）になるだけである。

          ∝ (1/√(τ^2+si^2))exp{-(<θi> - θ)^2/2(τ^2 + si^2)}

i について積を採れば、θ と τ を推定する周辺尤度の筈であるが、その対数を採っても式(9.15) の対数尤度とは異なり、最後の項が出てこない。
僕の理解はどこか間違っているのだろうと思う。

          ∝ ∑i {-(<θi> - θ)^2/(τ^2+si^2) + log(τ^2+si^2)} + log{∑i 1/(τ^2+si^2)} (9.15)

いずれにしても、random effect  モデルにおいては、重み因子 wi が τ に依存することになり、その τ^2 は施設間偏差によって増加した分散の追加分である。

     wi(τ) = 1/(τ^2 + si^2)
となり、推定効果 θ* は

     θ* = ∑i wi(τ)<θi>/∑i wi(τ)
分散は 1/wi(τ)
となる。ともあれ、問題は τ をどうやって決めるかである。
最尤推定は θ と τ の２つについて行わねばならず、非線形であるから、一般的には反復計算が必要となるが、今日の計算機の環境では造作もない。教科書には式が載っているが、僕はまだ理解していない。理解したら改定する。

・・・DerSimonian and Laird モデル というのは、反復計算を必要としない近似法であり、標準的に用いられている方法のようである。これは上記の計算とは独立に導かれている。
均一性の検定で使われた Q1 を考える（ここが丹後先生の導出法である）。

    Q1 = ∑i wi(<θi> - θ~)^2 = ∑i wi(<θi> - θ)^2 - (∑i wi)(θ~ - θ)^2

Q1 は、θi が全て等しいという帰無仮説の元では K-1次の χ2 分布に従う統計量であり、観測値だけから計算できる。
右辺第１項 Q は θi＝θ であるという帰無仮説の元で計算できて、K次の χ2 分布に従い、
第２項 Q2 は θ＝0 であるという帰無仮説の元で計算出来て、１次の χ2 分布に従う。
ここで、DerSimomian 嬢は、<θi> は観測値であるが、確率変数として扱い、右辺を期待値（<θi> の分布で平均した値）に置き換えることにした。それはその分布における分散である。そもそも、<θi>の母集団平均 θ からの分散は si^2 + τ^2 であるというモデルであるから、第１項の期待値 E() は自明であり、

    E((<θi> - θ)^2) = si^2 + τ^2

となる。つまり、第１項の期待値は、（wi = 1/si^2 だから）∑i (1 + τ^2wi) である。
第２項の変形がなかなか判らなかったのだが、一日頭を冷やして、やっと判った。

  θ~ = ∑i wi<θi>/∑i wi を代入して、自乗を開けば、異なる施設 i 間の相関が無いというランダムモデルだから、

    E((θ~ - θ)^2) = E({∑i wi(<θi> - θ)/(∑i wi)}^2) = E(∑i wi^2(<θi> - θ)^2/(∑i wi)}^2)
                  = ∑i wi^2(si^2 + τ^2)^2/(∑i wi)}^2 = {∑i wi + τ^2∑i wi^2}/(∑i wi)}^2

と変形出来て、全体をまとめると、

    E(Q1) = ∑i (1 + τ^2wi) - 1 - τ^2∑i wi^2/(∑i wi) = K - 1 + τ^2{∑i wi - ∑i wi^2/(∑i wi)}

となる。これから τ^2 を決めることができる。

・・・要するに、観測された各施設での効果量 <θi> と、その期待分散 si^2 から、施設間での偏差が無いという帰無仮説に基づく χ2 統計量 Q1 が決まり、それで偏差があるかどうかが統計的に判断されるのであるが、偏差があると判断された場合には、その偏差 τ を取り込んだ random effect モデルを構成し、そのモデルから期待される統計量 E(Q1) と実際の統計量 Q1 が等しいという近似を置くことで、逆に偏差 τ を決めてしまうことができる、という筋書きになっている。なかなか巧妙である。