平均値差のメタ分析の演習
環境
sessioninfo::platform_info()## setting value
## version R version 3.5.2 (2018-12-20)
## os macOS Mojave 10.14.5
## system x86_64, darwin15.6.0
## ui X11
## language (EN)
## collate ja_JP.UTF-8
## ctype ja_JP.UTF-8
## tz Asia/Tokyo
## date 2019-06-28データ
- 作ったデータをCSV形式にして（xls, xlsx形式と比べて読み込みが早いため），Rで読み込む。
cl <- read.csv("Cooperative_Learning.csv")- 以下のようなデータを扱うこととなる
- 論文タイトル（1列目）は取り出さない
cl <- cl[c(2:9)]
cl## Author N1 N2 M1 M2 SD1 SD2 t.value
## 1 Gull (2015) 30 32 18.57 16.31 3.471 NA 13.980
## 2 Yapici (2016) 30 31 18.33 13.50 1.370 0.980 2.870
## 3 Tran (2014) 51 55 77.36 67.00 4.520 6.600 9.600
## 4 Opdecam and Everaert (2012) 117 236 11.11 9.94 NA NA -2.135
## 5 Lang (1987) 60 30 32.50 33.00 4.960 4.910 NA
## 6 Jalilifar (2010) 30 30 12.30 10.23 3.300 2.760 NA
## 7 ??? 100 50 134.39 125.61 21.840 30.050 NA
## 8 Bot and Eze (2016) 70 70 29.98 22.71 10.370 10.400 NA
## 9 Zarifi and Taghavi (2016) 25 25 60.56 57.20 5.450 3.200 NA
## 10 Khan and Akhtar (2017) 90 88 43.55 30.13 6.940 7.890 -12.040
## 11 Rabgay (2018) 90 88 NA NA NA NA NA
## 12 AgwuUdu (2017) 205 191 29.06 19.55 8.530 5.380 NA
## 13 Najmonnisa and Saad (2017) 65 63 74.45 55.10 10.540 12.755 NA研究ごとの効果量と分散を求める
平均値差が示されている研究
- 平均値差が示されている研究だけのデータを作る
dplyrパッケージが入っていない場合には，install.packages("dplyr")としてパッケージをインストールする
library(dplyr)cl.m <- dplyr::slice(cl, 2:3, 5:6, 8:10, 12:13)
cl.m <- dplyr::select(cl.m, 1:7)
cl.m## Author N1 N2 M1 M2 SD1 SD2
## 1 Yapici (2016) 30 31 18.33 13.50 1.37 0.980
## 2 Tran (2014) 51 55 77.36 67.00 4.52 6.600
## 3 Lang (1987) 60 30 32.50 33.00 4.96 4.910
## 4 Jalilifar (2010) 30 30 12.30 10.23 3.30 2.760
## 5 Bot and Eze (2016) 70 70 29.98 22.71 10.37 10.400
## 6 Zarifi and Taghavi (2016) 25 25 60.56 57.20 5.45 3.200
## 7 Khan and Akhtar (2017) 90 88 43.55 30.13 6.94 7.890
## 8 AgwuUdu (2017) 205 191 29.06 19.55 8.53 5.380
## 9 Najmonnisa and Saad (2017) 65 63 74.45 55.10 10.54 12.755効果量
研究によって用いられる尺度（テストなど）の性質が異なることから，複数研究間の平均値等を単純に比較することはできない。そのため，これらを比較可能とするように効果量を求める。イメージとしては，異なる2つの試験の結果を比較する際に偏差値を用いるのとほぼ同じである。
平均値差の効果量 \(( d )\) は，0.20で小さい効果，0.50で中程度の効果，0.80で大きい効果と解釈すると言われているが，研究の文脈によって解釈が異なるので注意すること。
効果量\(d\) は以下の式で求められる。
\[ d=\frac{M_1-M_2}{\sigma _{pooled}} \]
なお，\(\sigma _{pooled}\) は母集団標準偏差のことであり，2群間の標準偏差の平均のようなものであり，以下の式によって求められる。
\[ \sigma _{pooled}=\sqrt{\frac{(n_1-1) s_1^2+(n_2-1) s_2^2}{n_1+n_2-2}} \]
効果量の分散
統合の対象となる研究はいずれもケース数（対象者数）が異なる。この違いを適切に反映させて効果量を統合させる必要があるため，その重みを求める必要がある。この重みは分散の逆数で求められるため，効果量の分散を求めておく。効果量の分散\(( 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} \]
効果量とその分散をRで求める
上記のように効果量と分散を求めるとメタ分析が可能となるが，研究例が多い場合には，統計パッケージを用いて一度に求めた方が効率的である。例題データに対して2群の平均値差の効果量とその分散を求める際には，以下の手順で行う。
metaforパッケージを読み込むmetaforパッケージが入っていない場合には，install.packages("metafor")としてパッケージをインストールする
library(metafor)metaforを用いて，各研究の効果量とその分散を求める。escalcが効果量等を求めるための命令であり，SMDとは2群の平均値差，n1i，n2iはそれぞれ第1群，第2群の人数，m1i，m2iはそれぞれ第1群，第2群の平均通過率，sd1i，sd2iははそれぞれ第1群，第2群の標準偏差を指している。data=で効果量等を求めるデータの名前を指定している。append=TRUEという命令を加えることで，結果をデータフレームとして（次の分析に使えるように）cl.m.resという名前で格納させる。
cl.m.res <- escalc(measure = "SMD",
n1i = N1, n2i = N2,
m1i = M1, m2i = M2,
sd1i = SD1, sd2i = SD2,
data = cl.m, append = TRUE)
cl.m.res## Author N1 N2 M1 M2 SD1 SD2 yi
## 1 Yapici (2016) 30 31 18.33 13.50 1.37 0.980 4.0144
## 2 Tran (2014) 51 55 77.36 67.00 4.52 6.600 1.8058
## 3 Lang (1987) 60 30 32.50 33.00 4.96 4.910 -0.1003
## 4 Jalilifar (2010) 30 30 12.30 10.23 3.30 2.760 0.6716
## 5 Bot and Eze (2016) 70 70 29.98 22.71 10.37 10.400 0.6962
## 6 Zarifi and Taghavi (2016) 25 25 60.56 57.20 5.45 3.200 0.7400
## 7 Khan and Akhtar (2017) 90 88 43.55 30.13 6.94 7.890 1.7997
## 8 AgwuUdu (2017) 205 191 29.06 19.55 8.53 5.380 1.3210
## 9 Najmonnisa and Saad (2017) 65 63 74.45 55.10 10.54 12.755 1.6464
## vi
## 1 0.1977
## 2 0.0532
## 3 0.0501
## 4 0.0704
## 5 0.0303
## 6 0.0855
## 7 0.0316
## 8 0.0123
## 9 0.0418yiとviという列が追加されているが，これらの列の情報が平均効果量を求めるのに必要である。
平均値差が示されておらず，\(t\)値が示されている研究
- ここで用いているHedges の\(g\)とt値とは以下のような関係がある
\[ t = \sqrt{\frac{n_1 n_2}{n_1+n_2}}g \]
- したがって，\(t\)値からHedgesの\(g\)は以下の式で求められる
\[ g = \frac{t}{\sqrt{\frac{n_1 n_2}{n_1+n_2}}} \]
- ここで，
clというデータから，平均値差が示されておらず，t値が示されている研究を取り出す。
cl.t <- dplyr::slice(cl, 1,4)
cl.t## Author N1 N2 M1 M2 SD1 SD2 t.value
## 1 Gull (2015) 30 32 18.57 16.31 3.471 NA 13.980
## 2 Opdecam and Everaert (2012) 117 236 11.11 9.94 NA NA -2.135- 各研究の効果量と分散を求める。分散は上に示した式で求められる。
- Gull (2015)の場合は以下の通りである。
# 効果量
Gull.d <- cl.t[1,8] / sqrt((cl.t[1,2] * cl.t[1,3]) / (cl.t[1,2] + cl.t[1,3]))
Gull.d## [1] 3.552773# 分散
Gull.vd <- ((cl.t[1,2] + cl.t[1,3]) / (cl.t[1,2] * cl.t[1,3])) + (Gull.d^2 / (2*(cl.t[1,2] + cl.t[1,3])))
Gull.j <- 1 - ((3)/(4 * (cl.t[1,2] + cl.t[1,3] -2) - 1))
Gull.vg <- Gull.j^2 * Gull.vd
Gull.vg## [1] 0.1622246- Opdecam and Everaert (2012)の場合は以下の通りである。
# 効果量
Opd.d <- cl.t[2,8] / sqrt((cl.t[2,2] * cl.t[2,3])/(cl.t[2,2] + cl.t[2,3]))
Opd.d## [1] -0.2413995# 分散
Opd.vd <- ((cl.t[2,2] + cl.t[2,3]) / (cl.t[2,2] * cl.t[2,3])) + (Opd.d^2 / (2*(cl.t[2,2] + cl.t[2,3])))
Opd.j <- 1 - ((3)/(4 * (cl.t[2,2] + cl.t[2,3] -2) - 1))
Opd.vg <- Opd.j^2 * Opd.vd
Opd.vg## [1] 0.01281187- これらの結果をデータテーブルして分析できるようにする
# 研究名を取り出す
cl.t.study <- cl.t[c(1)]
# 効果量と分散を行列にまとめる
cl.t.dv <- data.frame(rbind(c(Gull.d, Gull.vg), c(Opd.d, Opd.vg)))
colnames(cl.t.dv) <- c(names(cl.m.res[8:9]))
# 研究名をつけた行列にする
cl.t.dv <- cbind(cl.t.study, cl.t.dv)
cl.t.dv## Author yi vi
## 1 Gull (2015) 3.5527725 0.16222464
## 2 Opdecam and Everaert (2012) -0.2413995 0.01281187メタ分析を行う
データの整形
cl.m.resとcl.t.dvを結合する
# cl.m.res をyiとviだけにする
cl.m.dv_ <- cl.m.res[c(1, 8:9)]
cl.m.dv <- data.frame(cl.m.dv_)
# cl.t.dv を結合する
cl.dv <- rbind(cl.m.dv, cl.t.dv)
cl.dv## Author yi vi
## 1 Yapici (2016) 4.0143732 0.19768314
## 2 Tran (2014) 1.8057897 0.05317115
## 3 Lang (1987) -0.1002764 0.05005586
## 4 Jalilifar (2010) 0.6716283 0.07042571
## 5 Bot and Eze (2016) 0.6962348 0.03030265
## 6 Zarifi and Taghavi (2016) 0.7400398 0.08547659
## 7 Khan and Akhtar (2017) 1.7997387 0.03157323
## 8 AgwuUdu (2017) 1.3209813 0.01231692
## 9 Najmonnisa and Saad (2017) 1.6464386 0.04184654
## 10 Gull (2015) 3.5527725 0.16222464
## 11 Opdecam and Everaert (2012) -0.2413995 0.01281187変量効果モデルによる統合
研究を統合し平均効果量を求める場合，固定効果モデルか変量効果モデルかを選択することとなる。教科書では固定効果モデルが用いられているが，教育研究の場合には変量効果モデルを用いた方がよいと考えられる。その理由は以下の通りである。
- 固定効果モデル
- 統合対象の各研究の効果の背後には共通の真値がある。
- 各研究の効果量の誤差は標本誤差(sampling error)だけによるものである。
- 変量効果モデル
- 統合対象の各研究の効果量は研究ごとに変動する（共通の真値なし）。
- 研究間分散は標本誤差だけでは説明できない（教育研究の場合には年齢や文脈で介入の効果が異なることが多い）。
変量モデルによる平均効果量を求める
上記の手続きによりcl.dvとして格納された各々の研究に対する効果量とその分散の情報を用いて，平均効果量などを求める。変量モデルによって平均効果量を求めるには，以下のようにrma関数を用いる。
res <- rma(yi, vi, data = cl.dv, method = "REML")
res##
## Random-Effects Model (k = 11; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 1.6432 (SE = 0.7642)
## tau (square root of estimated tau^2 value): 1.2819
## I^2 (total heterogeneity / total variability): 97.82%
## H^2 (total variability / sampling variability): 45.82
##
## Test for Heterogeneity:
## Q(df = 10) = 285.0669, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 1.4131 0.3942 3.5847 0.0003 0.6405 2.1858 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1メタ分析の結果をフォレストプロットで表現する
メタ分析の結果は フォレストプロットで表現するとわかりやすい。フォレストプロットを描くには，先に得た出力resを用いる。
forest(res)
しかし，これではわかりにくいので，以下のように出力の命令を工夫する。このコードを用いた出力は以下の通りとなる。
forest(res, slab = paste(cl.dv$Author),
xlim = c(-8, 10), cex = 1.00)
#op <- par(cex = 0.7, font = 2)
text(7.0, 13, "d", pos = 2)
text(9.8, 13, "[95% CI]", pos = 2)
text(-8.0, 13, "Author", pos = 4)
課題
以下の課題に取り組み，7月9日(火)の講義で発表すること。資料は受講院生数分用意すること。
- 上記の分析を，自身の環境で再現し，メタ分析の方法を理解すること。
- 以下の要件を満たす発表資料を作成すること。
- [問題と目的] 何についてのメタ分析なのかを説明する。
- [方法] 文献の収集基準(用いたデータベースと検索式)を明らかにする。
- [方法] 分析対象とした文献の選択基準を明らかにする。
- [結果] メタ分析の結果を示す。
- [考察] 結果から何が言えるのかを考察する。
- [感想] 自身の研究にどのように役立つかといったことをはじめとした，感想を述べる。