Chap6_Yanamoto

6.1効果量の考え方

効果量の報告が必要な理由

p値の算出にはサンプルサイズが大きく影響する。
サンプルサイズが大きくなればなるほど統計的な有意差が出てきやすくなる。

#サンプルサイズ，平均値，標準偏差を指定してデータを作成する関数
#gendat関数・・・乱数でデータ作成  n　サンプルサイズ　mu 平均値　sigma 標準偏差
# 指定した平均値 mu と標準偏差 sigma を持つ n 個の正規乱数を発生させる
gendat <- function(n, mu=0, sigma=1)
{
  x <- rnorm(n)　　 # n 個の正規乱数を発生
  return((x-mean(x))/sd(x)*sigma+mu)  # 基準化する
}
#データセットAの作成（サンプルサイズ，平均値，標準偏差の順）
a <- gendat(50, 30.00, 10.00)
b <- gendat(50, 32.00, 10.00)
#データセットBの作成（サンプルサイズ，平均値，標準偏差の順）
x<-gendat(500, 30.00, 10.00)
y<-gendat(500, 32.00, 10.00)

#t検定（データセットA）
t.test(a, b)

## 
##  Welch Two Sample t-test
## 
## data:  a and b
## t = -1, df = 98, p-value = 0.3198
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -5.968935  1.968935
## sample estimates:
## mean of x mean of y 
##        30        32

#t検定（データセットB）
t.test(x, y)

## 
##  Welch Two Sample t-test
## 
## data:  x and y
## t = -3.1623, df = 998, p-value = 0.001613
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.2410952 -0.7589048
## sample estimates:
## mean of x mean of y 
##        30        32

BのセットはAのデータセットのサンプルサイズを10倍にしたもの

#　箱ひげ図と蜂群図の描画のめにデータセットの変換score.A <- c(a, b)
score.A <- c(a, b)
group.A <- factor(c(rep("1(n=50)", length(a)), rep("2(n=50)", length(b))))
score.B <- c(x, y)
group.B <- factor(c(rep("1 (n = 500)", length(x)), rep("2 (n = 500)", length(y))))

#グラフを2つ並べて表示する設定
par(mfrow = c(1, 2))
#箱ひげ図の作図
library("beeswarm")
par(family = "HiraKakuProN-W3") #日本語のラベルの文字化けを防ぐ
boxplot(score.A ~ group.A, ylim = c(0, 70), main = "データセットA", xlab = "指導法", ylab = "score")
beeswarm(score.A ~ group.A, ylim = c(0, 70), pch = 16, cex = 0.5, add = TRUE)
par(family = "HiraKakuProN-W3")
boxplot(score.B ~ group.B, ylim = c(0, 70), main = "データセットB", xlab = "指導法", ylab = "score")
beeswarm(score.B ~ group.B, ylim = c(0, 70), pch = 16, cex = 0.5, add = TRUE)

2つの図を見ると平均値は違いがない。データセットBでのみ有意差が見られたのはサンプルサイズが大きいため
単純に「p値が小さいほど，差が大きい」と判断するのは誤った考え方
APAでは効果量の報告が推奨されている（2020年版）

6.2 効果量の大きさの基準

Cohen(1988) 効果量小d=0.2, 効果量中d=0.5, 効果量大d=0.8
η（イータ二乗）は、要因が小さい場合に使う。/ 偏イータ　partialη二乗は要因が多い場合に使う。

6.3 Rによる効果量の算出

#パッケージの読み込み
install.packages("compute.es",dependencies = TRUE, repos = "http://cran.rstudio.com/")

## 
## The downloaded binary packages are in
##  /var/folders/45/9b3fn7gn0fx8h4sv9xdpxdwm0000gn/T//Rtmp2CzZWz/downloaded_packages

library("compute.es")

#データセットAの効果量を算出
mes(mean(a), mean(b), sd(a), sd(b), n.1 = 50, n.2 = 50)

## Mean Differences ES: 
##  
##  d [ 95 %CI] = -0.2 [ -0.59 , 0.19 ] 
##   var(d) = 0.04 
##   p-value(d) = 0.32 
##   U3(d) = 42.07 % 
##   CLES(d) = 44.38 % 
##   Cliff's Delta = -0.11 
##  
##  g [ 95 %CI] = -0.2 [ -0.59 , 0.19 ] 
##   var(g) = 0.04 
##   p-value(g) = 0.32 
##   U3(g) = 42.13 % 
##   CLES(g) = 44.42 % 
##  
##  Correlation ES: 
##  
##  r [ 95 %CI] = -0.1 [ -0.29 , 0.1 ] 
##   var(r) = 0.01 
##   p-value(r) = 0.32 
##  
##  z [ 95 %CI] = -0.1 [ -0.3 , 0.1 ] 
##   var(z) = 0.01 
##   p-value(z) = 0.32 
##  
##  Odds Ratio ES: 
##  
##  OR [ 95 %CI] = 0.7 [ 0.34 , 1.42 ] 
##   p-value(OR) = 0.32 
##  
##  Log OR [ 95 %CI] = -0.36 [ -1.08 , 0.35 ] 
##   var(lOR) = 0.13 
##   p-value(Log OR) = 0.32 
##  
##  Other: 
##  
##  NNT = -19.53 
##  Total N = 100

#データセットBの効果量を算出
mes(mean(x), mean(y), sd(x), sd(y), n.1 = 500,n.2 = 500)

## Mean Differences ES: 
##  
##  d [ 95 %CI] = -0.2 [ -0.32 , -0.08 ] 
##   var(d) = 0 
##   p-value(d) = 0 
##   U3(d) = 42.07 % 
##   CLES(d) = 44.38 % 
##   Cliff's Delta = -0.11 
##  
##  g [ 95 %CI] = -0.2 [ -0.32 , -0.08 ] 
##   var(g) = 0 
##   p-value(g) = 0 
##   U3(g) = 42.08 % 
##   CLES(g) = 44.38 % 
##  
##  Correlation ES: 
##  
##  r [ 95 %CI] = -0.1 [ -0.16 , -0.04 ] 
##   var(r) = 0 
##   p-value(r) = 0 
##  
##  z [ 95 %CI] = -0.1 [ -0.16 , -0.04 ] 
##   var(z) = 0 
##   p-value(z) = 0 
##  
##  Odds Ratio ES: 
##  
##  OR [ 95 %CI] = 0.7 [ 0.56 , 0.87 ] 
##   p-value(OR) = 0 
##  
##  Log OR [ 95 %CI] = -0.36 [ -0.59 , -0.14 ] 
##   var(lOR) = 0.01 
##   p-value(Log OR) = 0 
##  
##  Other: 
##  
##  NNT = -19.53 
##  Total N = 1000

mes関数では2群間の平均値，標準偏差，サンプルサイズがあれば効果量を算出できる
［95 %CI］は効果量が正しい場合に，同じ実験を繰り返したら95回はこの範囲に値が収まるという意味
データセットAの信頼区間とデータセットBの信頼区間の比較してみましょう ※ 信頼区間が０をまたいでいる場合，２つのグループに差がないと解釈する

6.4 ４章のデータから効果量を計算

#4章の対応のないt検定のデータから効果量を算出

mes(60.24, 72.27, 15.81, 16.11, 33, 37)

## Mean Differences ES: 
##  
##  d [ 95 %CI] = -0.75 [ -1.24 , -0.27 ] 
##   var(d) = 0.06 
##   p-value(d) = 0 
##   U3(d) = 22.56 % 
##   CLES(d) = 29.71 % 
##   Cliff's Delta = -0.41 
##  
##  g [ 95 %CI] = -0.74 [ -1.23 , -0.26 ] 
##   var(g) = 0.06 
##   p-value(g) = 0 
##   U3(g) = 22.81 % 
##   CLES(g) = 29.92 % 
##  
##  Correlation ES: 
##  
##  r [ 95 %CI] = -0.36 [ -0.55 , -0.13 ] 
##   var(r) = 0.01 
##   p-value(r) = 0 
##  
##  z [ 95 %CI] = -0.37 [ -0.61 , -0.13 ] 
##   var(z) = 0.01 
##   p-value(z) = 0 
##  
##  Odds Ratio ES: 
##  
##  OR [ 95 %CI] = 0.26 [ 0.11 , 0.62 ] 
##   p-value(OR) = 0 
##  
##  Log OR [ 95 %CI] = -1.37 [ -2.25 , -0.49 ] 
##   var(lOR) = 0.2 
##   p-value(Log OR) = 0 
##  
##  Other: 
##  
##  NNT = -6.91 
##  Total N = 70

有意差ありの結果でも効果量の結果から大きな差があることが確認できる

#4章の対応のあるt検定のデータから効果量を算出
#対応なしの場合のdの計算
res <- mes(67.33, 74.43, 9.66, 8.98, 30, 30)

## Mean Differences ES: 
##  
##  d [ 95 %CI] = -0.76 [ -1.29 , -0.24 ] 
##   var(d) = 0.07 
##   p-value(d) = 0.01 
##   U3(d) = 22.32 % 
##   CLES(d) = 29.52 % 
##   Cliff's Delta = -0.41 
##  
##  g [ 95 %CI] = -0.75 [ -1.27 , -0.23 ] 
##   var(g) = 0.07 
##   p-value(g) = 0.01 
##   U3(g) = 22.62 % 
##   CLES(g) = 29.76 % 
##  
##  Correlation ES: 
##  
##  r [ 95 %CI] = -0.36 [ -0.56 , -0.12 ] 
##   var(r) = 0.01 
##   p-value(r) = 0.01 
##  
##  z [ 95 %CI] = -0.38 [ -0.64 , -0.12 ] 
##   var(z) = 0.02 
##   p-value(z) = 0.01 
##  
##  Odds Ratio ES: 
##  
##  OR [ 95 %CI] = 0.25 [ 0.1 , 0.65 ] 
##   p-value(OR) = 0.01 
##  
##  Log OR [ 95 %CI] = -1.38 [ -2.33 , -0.43 ] 
##   var(lOR) = 0.24 
##   p-value(Log OR) = 0.01 
##  
##  Other: 
##  
##  NNT = -6.87 
##  Total N = 60

res

##   N.total n.1 n.2     d var.d   l.d   u.d  U3.d  cl.d cliffs.d pval.d     g
## 1      60  30  30 -0.76  0.07 -1.29 -0.24 22.32 29.52    -0.41   0.01 -0.75
##   var.g   l.g   u.g  U3.g  cl.g pval.g     r var.r   l.r   u.r pval.r fisher.z
## 1  0.07 -1.27 -0.23 22.62 29.76   0.01 -0.36  0.01 -0.56 -0.12   0.01    -0.38
##   var.z   l.z   u.z   OR l.or u.or pval.or   lOR l.lor u.lor pval.lor   NNT
## 1  0.02 -0.64 -0.12 0.25  0.1 0.65    0.01 -1.38 -2.33 -0.43     0.01 -6.87

この場合は，式(6.2）を使用
res[, “d”]は効果量dを示す

d.val <- res[, "d"] / sqrt(2*(1-0.6487101))
res[, 4]

## [1] -0.76

d.val

## [1] -0.9067045

library(psych)
d.ci(d.val, n1=30)

##          lower     effect      upper
## [1,] -1.327785 -0.9067045 -0.4745329

まとめ

効果量はp値のように実質的な差や効果を誤って解釈する可能性を避けるためのもの
かつ記述統計（平均値，標準偏差）に近い形で検討できる指標
結果を読み手に正確にわかりやすく伝えるために効果量を記述しよう

おまけ

検定力(power)
第一種の誤り（α）・・・本当は有意差がないのに有意差があると判断する確率
第二種の誤り（β）・・・実際には有意差があるのに有意差はないと判断する確率
検定力の算出（1-β)
RでGpowerの算出(Rmdで読み込めなかったので直接お伝えします。)

Chap6_Yanamoto

Daichi Yanamoto

6/30/2022

6.1効果量の考え方

効果量の報告が必要な理由

6.2 効果量の大きさの基準

6.3 Rによる効果量の算出

6.4 ４章のデータから効果量を計算

まとめ

おまけ