母数モデルに続いて変量モデル（変量効果モデル・ランダム効果モデル）と言われるメタアナリシスの方法について、統計数理的な観点から見ていこうと思います。
大まかな流れは2021年の統計検定1級応用、医薬生物学の問4をみるとわかりやすいので、一度解いてみると良いのかもしれません。
解説はこちらの記事に記載しました。
【統計応用・医薬生物学】2021年統計応用(医薬生物学)の解答例 問4【統計検定1級対策】
では概要を順番に見ていきましょう。
変量モデルの概要
①変量モデルの式
1~K個の研究が存在し、統合すると仮定します。
以前の記事では母数モデルの式を記載しましたが、その際モデルにおける真の効果量をθとしてそれが互いに等しいとしました。
変量モデルでは各研究における効果量\(\theta_i\)が研究ごとに異なる確率変数であると仮定します。
②漸近的正規分布となるように適切な関数で変換
母数モデルの場合と同様に、効果量θを適切な関数による変換することで、正規分布に従うようにします。分散\(\tau^2\)として以下のようになります。
\(f(\theta_i)\sim N(f(\theta), \tau^2)\)
個々の研究の効果量は分散\(\tau^2\)でばらつくことを想定します。さらに個々の推定値は誤差をもって測定されるため、母数モデルと同様に分散\(\s_i^2)となる正規分布に漸近的に従うことを仮定して
\(f(\hat\theta_i)\sim N(f(\theta), s_i^2 \tau^2)\)
となります。
③制限付き最尤推定などでパラメータの推定を行う
ここからパラメータの最尤推定を行うわけですが、母数モデルと異なり、固定効果(fixed effects)と変量効果(random effects)のパラメータが存在するため、推定が容易ではありません。
そこで、制限付き最尤推定(REstricted Maximum Likelihood;REML)を行うのが良いとされています。
基本的な考えと流れは前回の記事に記載しました。
【統計応用・医薬生物学】制限付き最尤推定量(REML; restricted maximum likelihood)の導出【統計検定1級対策】
制限付き最尤推定はどうしても反復計算を必要とするため、モーメント法を利用して推定を行うDersimonian-Lairdの方法もあります。詳細は後述していきます。
④重み付けをして、関数変換した効果量の推定値を統合
ここからは母数モデルと同様です。変量モデルでも同様に分散の逆数を重みとして推定値の統合を行います。つまり
\(w_i=\frac{1}{s^2_i＋\tau^2}\)
として統合した効果量は③で計算された推定値を上記のτに代入して
\(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％信頼区間は
\(f(\hat\theta)±1.96\sqrt{\frac{1}{\sum w_i}}\)
を同様に逆関数で戻すことで計算できます。
変量モデルの具体的な方法
さて、では変量モデルの具体的な方法を見てみます。
以下では毎回f(θ)と書くと煩雑なので各研究の結果を\(\hat\theta_i(i=1,2,…,K)\)として
\(\hat\theta_i\sim N(\theta, s_i^2＋\tau^2)\)
とします。
DerSimonian-Lairdの方法
まず反復計算のいらないモーメント法を用いた推定量の算出方法です。ここでは母数モデルの記事で紹介した均質性の検定統計量(Qと書かれることが多いと思います)を用いて、モーメント法を使っていきます。
均質性の検定についてはこちら
【統計応用・医薬生物学】メタアナリシスにおける母数モデルの統計数理と代表的な推定量【統計検定1級対策】
均質性の検定は帰無仮説として全ての研究の\(\theta_i\)が同様であることを仮定しているため、ここでは母数モデルと同様に重みを\(w_i=\frac{1}{s_i}\)とします。効果量θの最尤推定値を使って検定統計量Qを表すと
\(Q=\sum w_i(\theta_i-\hat\theta_{MLE})^2\)
となります。続いて、この期待値を求めてモーメント法を適用します。この期待値の計算では以下の少しテクニカルな展開をします。REMLの導出でも用いたお馴染みの式変形です。
\(\sum w_i(\theta_i-\hat\theta_{MLE})^2=\sum w_i(\hat\theta_i-\theta)\\=\sum w_i(\hat\theta_i-\theta)^2-\sum w_i(\hat\theta_{MLE}-\theta)^2\)
展開して和の記号も展開していけば等式が成り立つことがわかると思います。これを使いますと、均質性の検定統計量Qの期待値は
\(E[Q]=\sum w_iE[(\hat\theta_i-\theta)^2]-\sum w_iE[(\hat\theta_{MLE}-\theta)^2]\\=\sum w_iV(\hat\theta_i)-(\sum w_i)V(\hat\theta_{MLE})\)
となります。それぞれの分散を考えますと
\(V(\hat\theta_i)=\frac{1}{w_i}＋\tau^2\)
\(V(\hat\theta_{MLE})=V(\frac{\sum w_i\hat\theta_i}{\sum w_i})\\=\frac{1}{(\sum w_i)^2}V(\sum w_i\hat\theta_i)\\=\frac{1}{(\sum w_i)^2}\sum \{w_i^2(\frac{1}{w_i}＋\tau^2)\}\\=\frac{1}{\sum w_i} \frac{\sum w_i^2\tau^2}{(\sum w_i)^2}\)
と表されます。よって、Qの期待値は
\(E[Q]=\sum w_i(\frac{1}{w_i}＋\tau^2)-(\sum w_i)\{\frac{1}{\sum w_i}＋\frac{\sum w_i^2\tau^2}{(\sum w_i)^2}\}\\=K-1＋\tau^2(\sum w_i-\frac{\sum w_i^2}{\sum w_i})\)
となります。あとはモーメント法を使えば
\(\tau^2=max(0, \frac{Q-(K-1)}{\sum w_i-\frac{\sum w_i^2}{\sum w_i}})\)
となります。なお、\(\tau^2\)は負の値を取らないため、\(Q_1-(K-1)\lt0\)の時には\(tau^2\)=0となります。この時、個々の研究結果の真の効果量のばらつきが0ということになるので、当然ながら母数モデルと一致することになります。
REMLを求める方法
Desimonian-Lairdの方法はモーメント法で求めていましたが、制限付き最尤推定を使って推定する方法もあります。
今回はパラメータにθとτがありますが、まずθを局外パラメータとしてτを推定していきます。まず、周辺尤度を計算するため、尤度関数をパラメータθをくくり出せるように分解します。
\(L(\theta, \tau^2)=\Pi_{i=1}^K\{2\pi(s_i＋\tau^2)\}^{-\frac{1}{2}}exp\{-\sum\frac{(\hat\theta_i-\theta)^2}{2(s_i^2＋\tau^2)}\}\\=\Pi_{i=1}^K\{2\pi(s_i^2＋\tau^2)\}^{-\frac{1}{2}}exp\{-\sum\frac{(\hat\theta_i-\hat\theta_{MLE})^2＋(\hat\theta_{MLE}-\theta)^2}{2(s_i^2＋\tau^2)}\}\\=\Pi_{i=1}^K\{2\pi(s_i＋\tau^2)\}^{-\frac{1}{2}}exp\{-\sum\frac{(\hat\theta_i-\hat\theta_{MLE})^2}{2(s_i^2＋\tau^2)}\}×\sum\{2\pi(s_i^2＋\tau^2)\}^{-\frac{1}{2}}exp\{-\sum\frac{(\hat\theta_{MLE}-\theta)^2}{2(s_i^2＋\tau^2)}\}×\{2\pi\sum(s_i^2＋\tau^2)\}^{\frac{1}{2}}\)
となります。ここで以下の数式
\(\sum\{2\pi(s_i^2＋\tau^2)\}^{-\frac{1}{2}}exp\{-\sum\frac{(\hat\theta_{MLE}-\theta)^2}{2(s_i^2＋\tau^2)}\}\)
これは正規分布の再生性から
\(N(0,\sum(s_i^2＋\tau^2))\)
に従うため、周辺尤度をとるために積分すると1になることがわかります。よって
\(\int L(\theta, \tau^2)d\theta=\Pi_{i=1}^K\{2\pi(s_i＋\tau^2)\}^{-\frac{1}{2}}exp\{-\sum\frac{(\hat\theta_i-\hat\theta_{MLE})^2}{2(s_i^2＋\tau^2)}\}×\sum\{2\pi\sum(s_i^2＋\tau^2)\}^{\frac{1}{2}}\)
となります。
あとはここから対数尤度を求めて
\(logL(\tau^2)\propto-\sum log(s_i^2＋\tau^2)-\sum\frac{(\theta_i-\hat\theta_{MLE})^2}{s_i^2＋\tau^2}＋log\sum(s_i^2+\tau^2)\)
偏微分を行って
\(\frac{\partial}{\partial\tau^2}logL=-\sum\frac{1}{s_i+\tau^2}＋\sum\frac{(\hat\theta_i-\hat\theta_{MLE})^2}{(s_i^2＋\tau^2)^2}＋\frac{1}{\sum(s_i＋\tau^2)}\)=0
ここで\(\frac{1}{s_i^2＋\tau^2}=w_i\)として式を整理すると
\(\frac{\sum w_i^2}{\sum w_i}=\sum\{w_i+w_i^2(\hat\theta_i-\hat\theta_{MLE})^2\}\)
となります。あとはこの式に適当な初期値を入れて反復計算でしていくことでREMLの推定値を求めていきます。
参考文献：
コメントを残す