メタ分析

以下の内容ができるようになることがこの講義の重要な目的のひとつであるので，しっかりと習得すること。

例題データ

• このデータは地方自治体独自で実施された学力検査のうち，家庭学習に取り組む時間と学力検査得点との関係を検討したものの結果（公開されている情報）を用いて作られている。
• 公開されている情報の多くはグラフ程度であったため，なんとなく（物差しをあてがって推定）求めた平均点（M）と，経験則に基づきながら適当に乱数を発生させて標準偏差（SD）を生成して，メタ分析の実習ができるように作られている。
• 家庭学習に取り組む時間が一日1時間未満（hw_lo）と，1時間以上2時間未満(hw_hi)の児童生徒だけをデータ化している。

データ

コードブック

メタ分析の準備

メタ分析を行う際には

• ケース数（人数）
• 平均値
• 平均値の標準偏差

を用いて各々の研究についての 効果量 と 効果量の分散 を求める必要がある。

効果量

研究によって用いられる尺度（テストなど）の性質が異なることから，複数研究間の平均値等を単純に比較することはできない。そのため，これらを比較可能とするように効果量を求める。イメージとしては，異なる2つの試験の結果を比較する際に偏差値を用いるのとほぼ同じである。

平均値差の効果量 \( d \) は，0.20で小さい効果，0.50で中程度の効果，0.80で大きい効果と解釈すると言われているが，研究の文脈によって解釈が異なるので注意すること。

効果量d は以下の式で求められる。

\[ d=\frac{M_1-M_2}{\sigma _{pooled}} \]

なお， \( \sigma _{pooled} \) は母集団標準偏差のことであり，２群間の標準偏差の平均のようなものであり，以下の式によって求められる。

\[ \sigma _{pooled}=\sqrt{\frac{(n_1-1) s_1^2+(n_2-1) s_2^2}{n_1+n_2-2}} \]

たとえば，例題データの「柴絵山市」の「国語A」について，家庭学習時間の多少と学力テストの平均通過率との関係を検討するために効果量を求める場合には，以下の通りとなる。まず，プールされた標準偏差\( \sigma _{pooled} \)は，

\[ \sqrt{\frac{646\times 0.05^2+987\times 0.11^2}{646+987-2}}=0.09 \]

となり，効果量\( d \)は \[ d=\frac{0.871-0.886}{0.09}=-0.166 \] となる。

効果量の分散

統合の対象となる研究はいずれもケース数（対象者数）が異なる。この違いを適切に反映させて効果量を統合させる必要があるため，その重みを求める必要がある。この重みは分散の逆数で求められるため，効果量の分散を求めておく。効果量の分散\( V_g \)は

\[ V_g=J^2\times V_d \]

であり，\( V_d \)は

\[ V_d=\frac{n_1+n_2}{n_1 n_2}+\frac{d^2}{2(n_1+n_2)} \]

\( J \)は

\[ J=1-\frac{3}{4(n_1+n_2-2)-1} \]

である。したがって，先に求めた効果量\( d \)の分散を求めるには，まず\( V_d \)を以下の通り求め，

\[ V_d=\frac{646+987}{646\times 987}+\frac{-0.166^2}{2(646+987)}=0.0026 \]

\( J \)は

\[ J=1-\frac{3}{4(646+987-2)-1}=0.999 \]

となり，\( V_g \)は

\[ V_g=0.99^2\times 0.0026=0.0026 \] となる。

効果量とその分散をRで求める

上記のように効果量と分散を求めるとメタ分析が可能となるが，研究例が多い場合には，統計パッケージを用いて一度に求めた方が効率的である。例題データに対して2群の平均値差の効果量とその分散を求める際には，以下の手順で行う。

• 例題データを読み込む。なお，例題データは https://dl.dropboxusercontent.com/u/7515903/meta.csv にある。

meta_hw <- read.csv("meta.csv")

次に，metafor というパッケージを用いて，meta_hw というデータの効果量とその分散を求める。 escalc が効果量等を求めるための命令であり，SMD とは２群の平均値差， n1i，n2i はそれぞれ第1群，第2群の人数， m1i，m2i はそれぞれ第1群，第2群の平均通過率， sd1i，sd2i ははそれぞれ第1群，第2群の標準偏差を指している。 data= で効果量等を求めるデータの名前を指定している。 append=TRUE という命令を加えることで，結果をデータフレームとして（次の分析に使えるように） datという名前で格納させる。

library(metafor)

## Loading required package: Formula
## 
## Loading 'metafor' package (version 1.8-0). For an overview 
## and introduction to the package please type: help(metafor).

dat <- escalc(measure = "SMD", n1i = N_hw_lo, n2i = N_hw_hi, m1i = M_hw_lo, 
    m2i = M_hw_hi, sd1i = SD_hw_lo, sd2i = SD_hw_hi, data = meta_hw, append = TRUE)

なお，ここでできた dat というデータは以下のようになっている。 yi と vi という列が追加されているが，これらの列の情報が統合効果量を求めるのに必要である。

head(dat)

##   ID   Region Year Grade Subject N_hw_lo M_hw_lo SD_hw_lo N_hw_hi M_hw_hi
## 1  1 柴絵山市   23     6  Koku_A     646   0.871     0.12     987   0.886
## 2  2 柴絵山市   23     6  Koku_B     646   0.664     0.14     987   0.679
## 3  3 柴絵山市   23     6   San_A     646   0.643     0.05     987   0.671
## 4  4 柴絵山市   23     6   San_B     646   0.543     0.12     987   0.550
## 5  5 柴絵山市   23     9  Koku_A     269   0.683     0.05     721   0.700
## 6  6 柴絵山市   23     9  Koku_B     269   0.517     0.15     721   0.550
##   SD_hw_hi      yi     vi
## 1     0.10 -0.1384 0.0026
## 2     0.05 -0.1558 0.0026
## 3     0.11 -0.3071 0.0026
## 4     0.15 -0.0504 0.0026
## 5     0.13 -0.1490 0.0051
## 6     0.05 -0.3704 0.0052

メタ分析の実施

変量モデルによる平均効果量を求める

上記の手続きにより dat として格納された各々の研究に対する効果量とその分散の情報を用いて，平均効果量などを求める。変量モデルによって平均効果量を求めるには，以下のように rma 関数を用いる。

res <- rma(yi, vi, data = dat)

その結果を出力させると下記の通りとなり， estimate として出力されている部分を見ると，平均効果量は -0.2753 であることが分かる。ただし，この場合

\[ 第1群(家庭学習時間1時間未満)-第2群(家庭学習時間1時間以上) \]

として効果量を求めているため，推定された効果量の符号が負になっている。

res

## 
## Random-Effects Model (k = 43; tau^2 estimator: REML)
## 
## tau^2 (estimated amount of total heterogeneity): 0.0247 (SE = 0.0056)
## tau (square root of estimated tau^2 value):      0.1573
## I^2 (total heterogeneity / total variability):   99.13%
## H^2 (total variability / sampling variability):  115.00
## 
## Test for Heterogeneity: 
## Q(df = 42) = 4265.9061, p-val < .0001
## 
## Model Results:
## 
## estimate       se     zval     pval    ci.lb    ci.ub          
##  -0.2753   0.0245 -11.2360   <.0001  -0.3233  -0.2273      *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

メタ分析の結果をフォレストプロットで表現する

メタ分析の結果は フォレストプロット で表現するとわかりやすい。フォレストプロットを描くには，先に得た出力 res を用いる。

forest(res)

しかし，これではわかりにくいので，以下のように出力の命令を工夫する。このコードを用いた出力は以下の通りとなる。

forest(res, slab = paste(dat$Region, dat$Year, dat$Grade, dat$Subject, sep = ", "), 
    xlim = c(-1.5, 1), cex = 0.75)
op <- par(cex = 0.7, font = 2)
text(1, 45, "効果量[95% 信頼区間]", pos = 2)
text(-1.5, 45, "地域・年度・学年・教科", pos = 4)

相関研究の結果を平均値差の研究群に統合する

2群の平均値差の効果量と相関は互いに換算可能であるため，これらを統合したメタ分析を行うことも可能である。

例えば「梅磯県」だけ，家庭学習の取り組み時間と教科の学力テストの得点との関係を上表のように相関で報告していた場合，以下のようにして効果量とその分散を求める。

相関係数を標準化された平均値差に変換する

標準化された平均値差を提示している研究結果に，相関係数を提示している研究結果を統合してメタ分析ができるようにするには，相関係数を標準化された平均値差に変換する手続きが必要である。この手続きを山田・井上(2012)に沿って説明すると，以下の通りとなる。

相関係数を効果量に変換するには

\[ d=\frac{2r}{\sqrt{1-r^2}} \]

とし，その効果量の分散を求めるには，相関係数の分散を

\[ V_r=\frac{(1-r^2)^2}{n-1} \]

で求めたうえで，

\[ V_d=\frac{4V_r}{(1-r^2)^3} \]

で求める。

これらの式を用いて，先ほど求めた相関係数を効果量に変換し，また分散を求めると以下の通りとなる。

• 効果量

\[ d=\frac{2\times 0.14}{\sqrt{1-(0.14)^2}}=\frac{0.28}{\sqrt{1-0.14^2}}=\frac{0.28}{0.99}=0.28 \]

• 分散 \[ V_d=\frac{4\frac{(1-r^2)^2}{n-1}}{(1-r^2)^3}=\frac{4\times 0.00008}{0.963^3}=\frac{0.00032}{0.893}=0.0004 \]

相関から求めた効果量を統合して分析する

データを用いて効果量 yi と分散 vi を求めたら，以下の表のように表形式でまとめる。その際， 効果量の正負に注意 する。この場合は先の分析結果にあわせて， 追加する効果量の数値を全て負 にしている。

このデータを，先ほど使ったデータにまとめる。先ほど使ったデータをcsv形式で書き出し，上表の結果をスプレッドシートで追加する編集を行う。 太字 が 追加した部分 である。

これを新たに適当な名前（dat_2.csv）で保存し，それを再びRに読み込ませる。

dat_2 <- read.csv("dat_2.csv")

このデータ（dat_2）に対するメタ分析は，先述した手順と同様である。

res_2 <- rma(yi, vi, data = dat_2)
res_2

## 
## Random-Effects Model (k = 47; tau^2 estimator: REML)
## 
## tau^2 (estimated amount of total heterogeneity): 0.0229 (SE = 0.0050)
## tau (square root of estimated tau^2 value):      0.1514
## I^2 (total heterogeneity / total variability):   99.03%
## H^2 (total variability / sampling variability):  103.07
## 
## Test for Heterogeneity: 
## Q(df = 46) = 4388.7228, p-val < .0001
## 
## Model Results:
## 
## estimate       se     zval     pval    ci.lb    ci.ub          
##  -0.2701   0.0226 -11.9722   <.0001  -0.3143  -0.2258      *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

この結果についてのフォレストプレットを以下のように描画すると，出力は下図の通りとなる。

forest(res_2, slab = paste(dat_2$Region, dat_2$Year, dat_2$Grade, dat_2$Subject, 
    sep = ", "), xlim = c(-1.5, 1), cex = 0.75)
op <- par(cex = 0.7, font = 2)
text(1, 50, "効果量[95% 信頼区間]", pos = 2)
text(-1.5, 50, "地域・年度・学年・教科", pos = 4)

引用文献

山田剛史・井上俊哉 (2012). メタ分析入門：心理・教育研究の系統的レビューのために 東京大学出版会