論文で報告されていない相関係数をサルベージ！（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(標本サイズ, 平均値1, 平均値2, 標準偏差1, 標準偏差2, 報告されているt値)

使い方！

たとえば，

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値と分散共分散行列を返してくれる。自動で図示もするよ。こんなかんじ。

f:id:kusanagik:20150406224940p:plain

たとえば，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-test

data: 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なので，いれてサルベージしてやる。