2021年の過去問でもメタアナリシスを題材にした問題が出題されていました。
こうした問題は誘導するような小問が多いと思うので、流れが分かっていれば得点しやすいように思います。そこで、今回は母数モデルでの代表的なメタアナリシスの方法とその統計数理的な理論の流れを確認してみます。
統計数理的な細かい背景はなかなかネット上や論文でも得られるものが少なく、書籍はなおさら情報がなく、わかる範囲のものをまとめた備忘録となってしまっていますので何か良い情報源やまとめがありましたらご教授いただければと思います。
母数モデルの概要
母数モデルでは、以下のような流れで統合した効果量を推定していきます。
①各研究ごとの真の効果量θは一定で、各研究で分散が一定な測定誤差によってのみばらつきが出ていると考える
②推定する効果量が漸近的に正規分布となるように適切な関数で変換をする
③分散に従って重み付けをして、関数変換した効果量の推定値を統合する
④逆関数で求めたい統合した推定値を算出する
それぞれ順番に数式を確認していきます。
①母数モデルの式
まず、真の効果量の値をθとして
全部でn個となる各研究の効果量をそれぞれ
\(\theta_1,\theta_2,…,\theta_n\)
とします。
ここで母数モデルでは
\(\theta_1=\theta_2=…=\theta_n\)
となるような均質性を仮定します。逆にいえばこの均質性が仮定できない場合であれば母数モデルを用いることは適切でないといえます。
この値からそれぞれ一定のバラつきがあることを仮定して、統合するのが母数モデルです。
効果量の値として使われるのはオッズ比、リスク比が代表的です。
②漸近的正規分布となるように適切な関数で変換
求めたい効果量はそのままでは正規分布にならない場合や、平均・分散の推定がしにくい場合があります。そこで、効果量を適切な関数によって変換して推定しやすくします。
つまり、ある研究の真の効果量を\(\theta_i\)、関数をf(θ)とすると
\(f(\hat\theta_i)\sim N(f(\theta_i),s_i^2)\)
となるような変換です。
後述していきますが、ここで「効果量及び分散の推定をどのように行うか」でそれぞれの方法が分かれます。
③重み付けをして、関数変換した効果量の推定値を統合
先ほどの正規分布の分散が①で説明したバラつきを示すわけですが、効果量を統合する際にばらつきが大きい（つまり精度が悪い）方が重みが小さくなるようにするとより望ましい結果が得られると考えられます。
効果量を統合する時の重み\(w_i\)は分散の逆数をとって
\(w_i^2=\frac{1}{s_i^2}\)
となります。
統合された推定量\((f(\hat\theta)\))はこの重み付き平均を用いて
\(f(\hat\theta)=\frac{\sum_{i=1}^nw_i f(\hat\theta_i)}{\sum_{i=1}^n w_i}\)
となります。
なお、これは漸近的最尤推定量になっています。尤度関数から考えて導出してみましょう。先ほど述べていたようにf(θ)は漸近的に正規分布に従うことを踏まえて尤度関数を考えると
\(L(\theta)=\Pi_{i=1}^n\frac{1}{\sqrt{2\pi s_i^2}}exp\{-\frac{(f(\hat\theta_i)-f(\theta))^2}{2s_i^2}\}\)
となります。パラメータに関係しない定数を無視して対数尤度関数を考えると
\(logL(\theta)\propto\sum\frac{(f(\hat\theta_i)-f(\theta))^2}{2s_i^2}\\=\sum w_i(f(\hat\theta_i)-f(\theta))^2\)
となるので、パラメータで微分し、＝0として整理すると
\(\frac{\partial}{\partial\theta}logL=-2\sum w_i(f(\hat\theta_i)-f(\hat\theta))=0\\f(\hat\theta)=\frac{\sum w_if(\hat\theta_i)}{\sum w_i)}\)
先ほどの重み付き平均が算出されます。
④逆関数で統合された推定量を求める
もともと知りたかったθを関数変換したものになっているので、知りたい効果量の推定値になるように、最後に逆関数で変換します。
\(\hat\theta=f^{-1}(\frac{\sum_{i=1}^nw_i f(\hat\theta_i)}{\sum_{i=1}^n w_i})\)
こうして知りたかった推定量を求めることができます。
また、95%信頼区間についても、統合された推定量の分散を考えれば導出することができます。
\(V(\frac{\sum_{i=1}^nw_i f(\hat\theta_i)}{\sum_{i=1}^nw_i})\\=\frac{\sum w_i^2V(f(\hat\theta_i))}{(\sum w_i)^2} (w_iは定数であるため)\\=\frac{\sum w_i^2・\frac{1}{w_i}}{(\sum w_i)^2}\\=\frac{1}{\sum w_i}\)
よって95%信頼区間は
\(f(\hat\theta)±1.96\sqrt{\frac{1}{\sum w_i}}\)
これを逆関数で変換したものとなります。
母数モデルの具体的な方法(オッズ比)
ここからはオッズ比を求めることを例にして、具体的な方法を見ていきます。以下では統合する研究のi番目の結果が下図のようなものであったと仮定します。
上述した流れによりますと、ここから効果量を推定し、それを適切な関数で変換することで正規分布に従う形にしていきます。
オッズ比については対数をとった対数オッズ比が正規分布に従うことが知られています。
まずはどのような正規分布に従うのかをみていきます。
上の例において表の中のa,cがそれぞれ二項分布に従うと考えてみます。アウトカムが起きる確率を介入群ではp,コントロール群ではqとして、互いに独立として
\(a\sim Bin(n, p), c\sim Bin(m,q)\)
となります。すると、このモデルにおけるオッズ比ωは
\(\omega=\frac{\frac{p}{1-p}}{\frac{q}{1-q}}\)
となります。ここで、二項分布であることをふまえるとpおよびqの最尤推定量は
\(\hat p=\frac{a}{n}, \hat q=\frac{c}{m}\)
となります。p,qの2次元のパラメータを用いた関数であるωはそれぞれの最尤推定量を代入することで
\(\hat\omega=\frac{\frac{\hat p}{1-\hat p}}{\frac{\hat q}{1-\hat q}}\\=\frac{ad}{bc}\)
として推定します。この推定値は平均をωとする1変量正規分布になることが知られています*1。
実は2019年の統計検定1級・統計応用（医薬生物学）の問3で多項分布において、それを示す問題が出ています。難しくて初見では全く解けませんでした。
では、この推定値\(\hat\omega\)の期待値と分散を計算してみましょう。
まず、期待値ですが
\(E[\hat\omega]=E[\frac{\hat p}{1-\hat p}]E[\frac{1-\hat q}{\hat q}] (p,qは独立のため)\)
となります。
ここで、a,bが二項分布に従うことから
\(E[\hat p]=E[\frac{a}{n}]\\=p\)
同様にして、\(\hat q\)の期待値もqとなります。デルタ法を利用して上記の式を一次まで近似すると
\(E[\frac{\hat p}{1-\hat p}]E[\frac{1-\hat q}{\hat q}]\\=\frac{p}{1-p}\frac{q}{1-q}\\=\omega\)
となります。デルタ法の近似についてはこの後でも使われますのでわからなかったらこちらをご参照ください。
分散は対数オッズ比の分散をデルタ法を用いて計算し、そこから逆に計算すると簡単です。
\(V(log\hat\omega)=V(log\frac{\hat p}{1-\hat p})＋V(log\frac{\hat q}{1-\hat q})\)
ここで第1項をまず求めてみましょう。
\(V(\hat p)=\frac{p(1-p)}{n}\)に注意しつつデルタ法を用いて
\(V(log\frac{\hat p}{1-\hat p})\approx\{\frac{1}{p(1-p)}\}^2\frac{p(1-p)}{n}\\=\frac{1}{np(1-p)}\\=\frac{1}{np}＋\frac{1}{n(1-p)}\\\approx\frac{1}{a}＋\frac{1}{b} (a bが十分大きいとしてpを\hat pで推定)\)
となります。同様にすると
\(V(log\frac{q}{1-q})\approx\frac{1}{c}＋\frac{1}{d}\)
と求まります。
以上から
\(V(log\hat\omega)=\frac{1}{a}＋\frac{1}{b}＋\frac{1}{c}＋\frac{1}{d}\)
という綺麗な形になります。
ここで再度デルタ法を考えますと
\(V(log\hat\omega)\approx\{\frac{1}{\omega^2}\}V(\hat\omega)\)
となるので、近似的には
\(V(\hat\omega)=\omega^2(\frac{1}{a}＋\frac{1}{b}＋\frac{1}{c}＋\frac{1}{d})\)
となります。
そんなわけでオッズ比の推定量\(\hat\omega\)は近似的に平均ωと上記の分散の正規分布に従うと言えます。対数オッズ比についてデルタ法を用いて考えれば漸近的に
\(logOR\sim N(log\frac{ad}{bc}, \frac{1}{a}＋\frac{1}{b}＋\frac{1}{c}＋\frac{1}{d})\)
となります。
この辺りの結果は個々の研究についての漸近的な分布の話でしたが、メタアナリシスで知りたいのは全ての研究を統合した推定値とその分散です。それをどのように推定していくかは具体的な方法によって異なります。
対数オッズ比の特性をふまえたうえで、具体的な方法をみていきます。
漸近分散法
まずは漸近分散法と呼ばれる方法です。この方法では各研究のサンプルサイズが大きい仮定のもとで、個々の対数オッズ比が上で示したような正規分布に近似的に従うことを利用して計算していきます。
各研究の対数オッズ比の推定量は、ある研究の結果を添字iを用いて表現すると上で示したように
\(log\hat\omega_i=log\frac{a_id_i}{b_ic_i}\)
分散は
\(V(log\hat\omega_i)=\frac{1}{a_i}＋\frac{1}{b_i}＋\frac{1}{c_i}＋\frac{1}{d_i}\)
となります。
この分散の平方根（SE）の逆数を重みwとし、重みつき平均を求め、逆関数で変換すれば（対数の逆関数なのでexpをつける）統合された推定量となります。
\(w_i=\frac{1}{\sqrt{V(log\hat\omega_i)}}\)
\(\hat{統合した\omega}=exp(\frac{\sum_{i=1}^nw_ilog\omega_i}{\sum w_i})\)
信頼区間は点推定量±1.96SEとして求めることができます。それぞれの研究結果が漸近的に正規分布に従うことが前提なので、サンプルサイズが大きいことが必要となります。
Mantel-Haenszel法
続いてはマンテル-ヘンツェル法です。
個々の研究(i番目の研究ということで全て添字iをつけます)のオッズ比を
\(\omega_i=\frac{a_id_i}{b_ic_i}\)
として推定し、そのオッズ比の標準誤差と重みを
\(SE_i=\frac{N_i}{b_ic_i}\)
\(w_i=\frac{b_ic_i}{N_i}\)
で推定します。
他の方法と同様に、標準誤差の逆数を重みとして、重み付き平均を出すことで統合オッズ比の推定値を算出します。
\(\hat{統合した\omega}=\frac{\sum\frac{a_id_i}{N_i}}{\sum\frac{b_ic_i}{N_i}}\)
なお、信頼区間についてはRobins-Breslow-Greenland varianceと呼ばれる統合したオッズ比の分散を用いて計算されますが、さすがに統計検定には出てこないと思われるので本記事では省略します。
で、こちらも統計検定には出されないとは思うのですが、気になるのがこの重み\(\frac{b_ic_i}{N_i}\)はどこからやってきたのかということです。
『メタアナリシス入門』*2によるとMantel-Haenszel推定量は数学的な導出というより天性のセンスによるもの、とされており、その理由を数学的に探ることは難しそうです。ただ、その性質が良いのかは事後的に議論がなされているようで、いくつか例をみてみます。
まず言えるのはこの推定量も統合した値も大数の法則が成り立つとき共通オッズ比と近似される点です*3。
大数の法則が成り立てば
\(\frac{a_i}{n_i}\approx p, \frac{b_i}{n_i}\approx 1-p, \frac{c_i}{m_i}\approx q, \frac{d_i}{m_i}\approx 1-q\)
となるので、共通オッズ比をωとすれば
\(\omega_i=\frac{a_id_i}{b_ic_i}\\\frac{\frac{a_id_i}{nm}}{\frac{b_ic_i}{nm}}\\\approx\frac{p(1-q)}{(1-p)q}\\=\omega\)
となります。重みをつけた場合でも、各研究のオッズ比が共通オッズ比に近似されるため(\(\omega_1=\omega_2=…=\omega\))以下のようになります。
\(\frac{\sum w_i\omega_i}{\sum w_i}=\frac{\omega\sum w_i}{\sum w_i}\\=\omega\)
というわけで統合されたオッズ比も各研究における大数の法則が成立すれば近似的に共通オッズ比となることがわかりました。
また大数の法則とともに帰無仮説ω＝1が成立する場合には、分散が先ほどの漸近分散と一致することもわかります*2,3。ω=1のとき、p=qとなるので
\(V=\frac{N_i}{b_ic_i}\\=\frac{\frac{N_i}{n_im_i}}{\frac{b_ic_i}{n_im_i}}\\\approx(\frac{1}{n_i}＋\frac{1}{m_i})\frac{1}{(1-p)p}\\=\omega^2(\frac{1}{np(1-p)}＋\frac{1}{mq(1-q)})\)
となるので、漸近分散と一致することがわかりました。
また別の説明として条件なしの最尤推定量（仮定が色々とつきます）とするものもあります*4。
やり方としては先ほどと同様に二項分布モデルを仮定します。するとi番目の研究の結果が上記の表のようになった場合、その尤度関数は添字をiとすると
\(L=p_i^{a_i}(1-p_i)^{b_i}q_i^{c_i}(1-q_i)^{d_i}\)
となります。これを全ての研究についての尤度関数とすると
\(\Pi_{i=1}^np_i^{a_i}(1-p_i)^{b_i}q_i^{c_i}(1-q_i)^{d_i}\)
となります。これを対数尤度関数として求めていくのですが、今回最尤推定したいパラメータはあくまで共通オッズ比であるωですので、p、qの一方を消去します。どちらでも良いのですがpに犠牲になってもらいましょう。
\(p_i=\frac{\omega q_i}{\omega q_i＋(1-q_i)}\)
と変換できますので、これを先ほどの対数尤度関数に代入してひたすら変形していきます。またωの推定に関連のないような局外パラメータはいずれにしても微分の際に消えるので省いていきます。
\(logL=\sum\{a_ilog\frac{\omega q_i}{\omega q_i (1-q_i)}＋b_i log\frac{q_i-1}{\omega q_i＋(1-q_i)}\}\\=\sum\{a_ilog\frac{\omega q_i}{\omega q_i＋(1-q_i)}＋b_i log(q_i-1)＋b_ilog\frac{1}{\omega q_i＋(1-q_i)}\}\\=\sum\{a_ilog\omega＋a_ilogq_i＋b_ilog(q_i-1)-(a_i＋b_i)log(\omega q_i＋(1-q_i))\}\)
さて、ここまで変形したらωで微分して
\(\frac{\partial}{\partial\omega}logL=\sum\{\frac{a_i}{\omega}-(a_i＋b_i)\frac{\omega q_i}{\omega q_i (1-q_i)}\}=0\\\sum\frac{\hat\omega q_ia_i a_i-q_ia_i-(a_i b_i)\hat\omega q_i}{\hat\omega q_i (1-q_i)}=0\)
ここで分母の値が研究毎に大きく変わらない（\(q_i\)があまり変わらない）と仮定すると
\(\sum\{\hat\omega q_ia_i＋a_i-q_ia_i-(a_i＋b_i)\omega q_i\}=0\\\hat\omega=\frac{\sum a_i(1-q_i)}{\sum q_i(a_ib_i-a_i)}\)
このとき\(q_i=\frac{d_i}{c_i＋d_i}\)として代入して求めると
\(\hat\omega=\frac{\sum a_i\frac{d_i}{c_i＋d_i}}{\sum \frac{d_i}{c_i＋d_i}b_i}\\=\frac{\sum\frac{a_id_i}{c_i＋d_i}}{\sum\frac{b_ic_i}{c_i＋d_i}}\)
ここで介入群とコントロール群が\(c_i＋d_i:a_i＋b_i=1:r\)と比が一定であるとするなら定数1＋rをかけると研究の参加者数の総数を\(N_i\)として
\((c_i＋d_i)(1＋r)=N_i\)
となるので、分子分母に1 rをかけて
\(\hat\omega=\frac{\sum\frac{a_id_i}{N_i}}{\sum\frac{b_ic_i}{N_i}}\)
というわけで「介入:コントロールの比が一定である」「コントロール群でのイベント発症率(q)があまり変わらない」という仮定のもとに条件なしの最尤推定量となることがわかりました。
Peto法
Petoの方法はランダム化比較試験の結果をメタアナリシスとして統合する際に用いられる方法です。
対数オッズ比の近似値を以下の数値で推定します*2,5。
\(log\hat\omega=\frac{O_i-E_i}{V_i}\)
ここでのOは観測度数（介入群のイベント数）、Eは超幾何分布の期待度数（治療効果がない場合に期待される介入群のイベント数）、Vはその分散となります。
よって、先ほどの2×2表の記号と同様にすると
\(log\hat\omega=\frac{O_i-E_i}{V_i}=\frac{a_i-\frac{n_i(a＋c)}{N_i}}{\frac{n_im_i(a＋c)(b＋d)}{N_i^2(N_i-1)}}\)
となります。
重みは上記の分散を用いて
\(w_i=V_i\)
となるので統合したオッズ比は、重み付けをして統合した対数オッズ比を逆関数で変換して戻せば求められます。
\(\hat{統合した\omega}=exp(\frac{\sum V_i\frac{O_i-E_i}{V_i}}{\sum V_i})\\=exp(\frac{\sum(O_i-E_i)}{\sum V_i})\)
ここでなぜO-Eの単純な形で対数オッズ比が表現できるかと言うことですが、方法としてはロジスティック回帰モデルで条件付き最尤推定を行い、βをフィッシャースコア法で初期値β＝0（つまりオッズ比が1）で一次近似した際に得られるようです*6,7。
フィッシャースコア法からの流れはわかるのですが、条件付き最尤推定の尤度関数がなぜそのような形になるのか、というあたりは元文献さかのぼって漁り続けましたが、どうにも綺麗に導出できませんでしたので、どなたかうまく説明できる方がいたら教えていただきたいところです。いつか解けたら追記したいと思います・・・。
以上、なんとかわかる範囲で資料を寄せ集めてメタ解析・母数モデルの統計数理を勉強しましたが、やはり応用分野はこれといってまとまった書籍がないので対策がしにくいですね。今後またうまくまとめられるような資料が見つかって学びが進んだらまたリライトしたいと思います。
おまけ：均質性の検定の導出
ちなみに均質性の検定、つまり推定した個々のオッズ比\(\omega_i\)と共通オッズ比ωが等しくなることを帰無仮説\(H_0\)とした検定は上で述べた対数オッズ比が正規分布に従うことを利用して行うことができます*2。推定値は個々の方法によって異なりますが基本的に以下のようになります。
個々の研究の推定オッズ比を\(\omega_i\)、共通オッズ比をωとすると帰無仮説は
\(H_0: \omega_1=\omega_2=…=\omega\)となる
このとき以下の検定統計量は自由度１のカイ二乗分布に従う（対数オッズ比は正規分布に従い、また重みは各推定対数オッズ比の分散の逆数となっているため）
\((log\omega_i-log\omega)^2w_i\sim\chi^2_1\)
よってカイ二乗分布の和の再生性からすべての研究の推定オッズ比の和を用いた以下の検定統計量Qは自由度k-1のカイ二乗分布に従う。
\(\sum_{i=1}^k(log\omega_i-log\omega)^2w_i\sim\chi^2_{k-1}\)
これを検定すればオッズ比が共通となるかどうかの確認をすることができます。棄却されてしまう場合はその核研究の結果を統合してよいものかどうか疑わしくなります。
上記の話では標準正規分布に従う確率変数を二乗したものがカイ二乗分布に従うこと、カイ二乗分布の和の再生性を用いていますが、以前に証明を記事にしていますので、気になったらこちらを参照してください。
標準正規分布とカイ2乗分布・ガンマ分布の関係について、整理と証明【統計検定1級対策】
それぞれの推定方法に従ってωを代入すればよいので、例えばマンテルヘンツェル法による共通オッズ比の推定値を\(\hat\omega_{MH}\)とすれば
\(\sum_{i=1}^k(log\omega_i-log\hat\omega_{MH})^2w_i\sim\chi^2_{k-1}\)
となります。Petoの方法の場合はもう少し変形することができて
\(\sum_{i=1}^k(log\omega_i-log\hat\omega)^2w_i\\=\sum\{\frac{O_i-E_i}{V_i}-\frac{\sum(O_i-E_i)}{\sum V_i}\}^2V_i\\=\sum\{\frac{(O_i-E_i)^2}{V_i}-2\frac{O_i-E_i}{V_i}\frac{\sum(O_i-E_i)}{\sum V_i}+(\frac{\sum(O_i-E_i)}{\sum V_i})^2\}V_i\\=\sum\frac{(O_i-E_i)^2}{V_i}-\frac{\{\sum(O_i-E_i)\}^2}{\sum V_i}\sim\chi^2_{k-1}\)
となります。
(2022/5/14追記)
参考文献：
*1
BIO 231 cd Statistical Inference Supplementary Course Notes
講義資料のようですが、p.53あたりからオッズ比が漸近的に正規分布に従うことが示されています。
*2『メタアナリシス入門』
*3Jupyter Notebook Viewer黒木玄先生が公開されている資料です。参考にさせていただきました。
Robins-Breslow-Greenland varianceについてはこちらがわかりやすいです。なお、Page2の左中央”by terms involving Ψ is:”以後の数式は和の記号の括弧が抜けていると思います。
*5Beta blockade during and after myocardial infarction: an overview of the randomized trials – PubMed
*6BIAS IN THE PETO ONE-STEP ESTIMATOR FOR THE COMMON ODDS RATIO
*7Statistical methods in cancer research. Volume I – The analysis of case-control studies – PubMed
コメントを残す