論文で報告されていない相関係数をサルベージ!(Rの関数つき)
背景
対応のあるデータの論文で,あるある。
平均値と標準偏差のペアが報告されている,よしよし。
対応のあるt検定の検定統計量が報告されている,ふむふむ。
でも相関係数や分散共分散行列が報告されていない!なんでだよー。。。
たとえば,標準化効果量は普通のd(標準化平均差)と,差得点に基づく標準化平均差では値が異なる。もちろん,相関係数がわかれば,このふたつは変換できるけど,相関係数が報告されていないとそれはむり。むり。メタ分析にも影響するよね。
もともと分散共分散行列,平均があればデータはだいたいシミュレートできる。
R,MASSパッケージのmvrnorm関数とか。でも,報告されていないとなー。
指導法の検証でも,そもそも共分散は教育的にも大きな意味を持つので,報告したほうがいいとおもうの。切実に。切実に。
でも,まあいいや。以下の条件なら,報告されていない相関係数がわかるから。
- 平均値のふたつが報告されている
- 標準偏差のふたつが報告されている
- 対応のあるt検定の検定統計量tがある
- (3の代わりにp値と標本サイズが報告されている)
関数つくっちゃえ!
対応ありのt検定にもちいる,「差得点に基づく標準偏差」は,完全データがなくても,二群の標準偏差と相関係数からももとまる。簡単にいえば,対応のあるt検定には相関係数が使われているので,そこからサルベージできるかもだ。
数式から計算もできるのだろうけど,ちょっとめんどくさいので,以下のようなアルゴリズムで答えをえちゃえ。
- 相関係数を0.000から,0.999まで用意する
- 論文で報告されている平均値と標準偏差と,用意したすべての相関係数をもちいて対応のあるt検定の検定統計量を計算する
- 論文で報告されているt値との差分をすべてもとめる
- 最小の差分に対応する相関係数が答え!
簡単だ♪尤度といえなくもない。
サルベージ!サルベージ!
これをRの関数でつくっちゃえ。
salvager<-function(n,m1,m2,sd1,sd2,t){
md<-abs(m1-m2) #平均差
v1<-sd1^2 #分散1
v2<-sd2^2 #分散2
sr<-seq(0,.999,.001) #相関係数の用意
sigmad<-sqrt(v1+v2-(2*sr*sd1*sd2)) #差得点の標準偏差
simt<-md/sigmad*sqrt(n) #対応あるt値を計算
diff<-t-simt #報告値との差分
ref<-data.frame(sr,abs(diff))
mindiff<-min(ref[,2]) #最小の差分をもとめる
rr<-ref[,1][ref[,2]==mindiff] #対応する相関係数
plot(ref,type="l",xlab="r", ylab="Difference") #図示
points(rr,mindiff,cex=2,pch=4,col=2,lwd=3)
varmatrix<-matrix(c(v1,rr*sd1*sd2,rr*sd1*sd2,v2),2,2) #分散共分散行列
list("r"=rr,"t"=t,"var"=varmatrix) #結果
}
使い方は,
使い方!
たとえば,
salvager(100,20,30,10,10,15.8114)
$r
[1] 0.8$t
[1] 15.8114$var
[,1] [,2]
[1,] 100 80
[2,] 80 100
こんな感じで,サルベージした相関係数と,報告されたt値と分散共分散行列を返してくれる。自動で図示もするよ。こんなかんじ。
たとえば,MASSのmvrnormで,
mvrnorm(100,mu=c(20,30),Sigma=matrix(c(100,30,30,100),2,2),empirical=T)->dat
とか作って,
t.test(dat[,1],dat[,2], paired=TRUE)
Paired t-testdata: dat[, 1] and dat[, 2]
t = -8.4515, df = 99, p-value = 2.563e-13
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-12.347757 -7.652243
sample estimates:
mean of the differences
-10
t値が8.4515なので,いれてサルベージしてやる。
salvager(100,20,30,10,10,8.4515)
$`Salvaged r`
[1] 0.3$t
[1] 8.4515$var
[,1] [,2]
[1,] 100 30
[2,] 30 100
おうおう!分散共分散行列が再現出来てるね。
念のため,
cor(dat)
[,1] [,2]
[1,] 1.0 0.3
[2,] 0.3 1.0
よしよし。
メタ分析にも役立つかもだし,レビューでも役立つかも。