【R言語】得られたデータから効果量を計算する（パラメトリック二群比較編）

※2020年5月4日に効果量について訂正&追記しています

こんにちは。プログラミング初心者のえいこです。

前回は、"pwr"というパッケージを使って（事前に実験を行うための）サンプル数を計算することをしてみました。

"pwr"はサクッと検定力やサンプルサイズを計算するためにはすごく便利なパッケージです。

でも、得られたデータが正規分布ではない場合の効果量は求められないのが欠点...

ということで今回は、得られたデータから効果量をしていきたいと思います。

効果量とは、

「サンプルサイズに影響されない効果の大きさを表す指標」

です。こちらの記事にもまとめています。

※下の記述は2020年5月4日に追記しています。

私なりの理解では「効果量とは検出できる差の大きさ」で、データの差を肉眼で見る（large）か、虫眼鏡で見ている（medium）か、顕微鏡で見ている（small）かの違いを示したものという感じでとらえています。

肉眼では区別がつかなかった差を、虫眼鏡で見れば検出できるかもしれないし、顕微鏡で見れば顕微鏡で見える範囲を超えている差になっているかもしれない...という感じです。

以前、有意水準、検定力を固定して効果量を変えてサンプルサイズを計算した時に、largeよりもmediumの方が必要なサンプル数が多いことがわかりました。

ということは、より精度が高い検出を行う（肉眼ではなくて虫眼鏡や顕微鏡で見る）ためにはサンプル数を増やさなければならないということです。

ただ、サンプル数を増やすと（pの値は小さくなって）第一種の過誤を起こしやすくなるので注意しましょうという感じ。

統計の専門家ではないので、正しくないかもしれません。

※2020年5月4日に追記
この肉眼・虫眼鏡・顕微鏡理論は効果量ではなくて、有意水準の概念ではないかと最近思うようになりました。

では、効果量とは？効果量には大きく分けてd族・r族があります。
d族は平均の差の大きさを、r族は相関の大きさを示しているようです。

と２種類あるようなのですが、肉眼・虫眼鏡・顕微鏡の倍率など関係なく、平均の差の大きさが1 cmなのか1 mmなのか1 μmを表す指標なようです。

もう少し考えを整理して、追記するかもしれません。

最近は効果量の掲載を求める科学系の論文が増えてきているので（心理学の分野ではかなり前から取り入れられてきているようですが）、きちんと効果量と有意水準を区別できるようになっておきたいものです。

パラメトリックな二群比較の場合の効果量の計算

パラメトリックな二群検定って一言で言えばt検定のことです。

t検定の時の効果量は一般的に"r"が使われるようなので、今回は、rが計算できるようになっておきましょう。計算式はこちらの記事でも紹介しています。

rは検定総計量のtと、自由度のdfを使って計算します。

この時に作ったデータを使ってt検定をした時の、出力結果はこちらです。

	Welch Two Sample t-test

data:  xtreat and placebo
t = 6.7815, df = 2.6212, p-value =
0.009927
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.9292014 2.8641319
sample estimates:
mean of x mean of y 
5.713333  3.816667 

この結果は、t値とdfが出力されているので効果量をこの結果を使って計算できそうです。

t検定の結果からtの値とdfの値を取ってくる。

p値を取り出せたので、tの値もdfの値も取り出せるのでは？ということでt検定の結果をviewで見てみようと思います。

ttest<-t.test(data$value ~ data$treatment)
view(ttest)

を実行すると...

こんな感じの画面が出てきます。"statistic"というところにt値が、"parameter"というところに自由度が計算されて入っています。

ということは、"$"を使って持ってくれば良いじゃん！簡単♪簡単♪

ということでやってみます。

t<-ttest$statistic
t
       t 
-6.781503 
df<- ttest$parameter
df
     df 
2.621168 

なんだかデータの形が２行一列になっています。データの形を見てみると..."named num ~"と書かれています。

つまり、値に名前がついてしまっているということ。このまま数式に落とし込むことができません。

ちなみに、p値の場合どうでしょう？

 p<-ttest$p.value
 p
[1] 0.009926739

数値データだけになっています。

”named”を除くためにはどうしたら良いのか検索したてみたところ、こちらのサイトが見つかりました。

書かれている回答の中で、「More generally, use」と書かれていたのが"unname()"だったので、今回はunname()で消してみたいと思います。

> t<-unname(ttest$statistic)
> t
[1] -6.781503
> df<- unname(ttest$parameter)
> df
[1] 2.621168

数値データだけになりました！

これを計算に使っていきます。

t.effect<- sqrt(t^2/(t^2+df))
t.effect
[1] 0.9726652

効果量は、0.97となりました。

次回はノンパラメトリック２群比較の効果量の計算をしてみたいと思います。

それでは、また！