はじめに
PythonにはPythonMetaというメタ分析を行うパッケージもあるのですが、勉強の意味も兼ねて、今回はpandasでメタアナリシスを実装してみました。
メタアナリシスのモデルには主に Fixed Effect Model （固定効果モデル） と Random Effects Model （変量効果モデル） の2種類があるので、それぞれ実装してみます。
サンプルデータ
今回は、以下のデータを使って実験してみます。
import pandas as pd
import numpy as np
data = pd.DataFrame({
"g": [0.12, 0.23, 0.34, 0.45, 0.42, 0.39, 0.49, 0.65, 0.76, 0.87],
"V": [0.01, 0.04, 0.03, 0.02, 0.01, 0.02, 0.03, 0.04, 0.02, 0.01]
})
ここで g が Hedges'g （＝効果量）、Vが効果量の分散を表しています。
固定効果モデル
固定効果モデルではシンプルに効果量の分散の逆数をその項目の重みとします。したがって求める平均効果量はそれぞれの項目の効果量に重みをかけたものを全体の重みで割ったものとして計算されます。計算式は以下の通りです。Pythonでも3行で書けtしまいます。
# 固定効果モデル
data['W'] = 1 / data.V
data['Wg'] = data['g'] * data['W']
result = data['Wg'].sum() / data['W'].sum()
result
>> 0.4776
変量効果モデル
変量効果モデルでは重みを計算する部分が少し複雑になります。研究間の効果量のばらつきを考慮するために、一旦効果量の平均をとったうえで、それに対する各項目の偏差を重みの計算に組み込みます。計算式は以下の通りです。
g_hat = data.g.mean()
Q = (data.W * (data.g - g_hat)**2).sum()
data['W2'] = data.W ** 2
C = data.W.sum() - (data.W2.sum()/data.W.sum())
d = len(data) - 1
# 研究間分散
if (Q-d) > 0:
V_between = (Q - d) / C
else:
V_between = 0
data['V_str'] = data.V + V_between
data['W_str'] = 1 / data.V_str
result = (data.g * data.W_str).sum() / data.W_str.sum()
result
その他
固定効果モデルでも変量効果モデルでも、計算された重みの和から95％信頼区間を計算することができます。
std_err = np.sqrt(1/data.W_str.sum())
lwr = result - 1.96 * std_err
upr = result + 1.96 * std_err
[lwr, upr]