1 効果量

効果量1
効果量2
d族の効果量には「Cohenのd」と「Hedgesのg」があります。これらは2つの標本間の平均値の差を標準偏差で割って標準化したもので、2標本の平均値の差を比較する際に使います。すなわち、2標本の平均値がどれだけ離れているかを表します。

1.1 データ作成関数

gendat <- function(n, mu = 0, sigma = 1)
{
  x <- rnorm(n)
  return((x - mean(x)) / sd(x) * sigma + mu)
}

1.1.1 データセットAの作成

a <- gendat(50, 30.00, 10.00)
b <- gendat(50, 32.00, 10.00)

1.1.2 データセットBの作成

x <- gendat(500, 30.00, 10.00)
y <- gendat(500, 32.00, 10.00)

1.1.3 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

#データセット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

1.2 因子型に変換

score.A <- c(a, b)
group.A <- factor(c(rep("1 (n = 50)", length(a)), rep("2 (n = 50)", length(b)))) 
str(group.A);str(score.A)

##  Factor w/ 2 levels "1 (n = 50)","2 (n = 50)": 1 1 1 1 1 1 1 1 1 1 ...

##  num [1:100] 24.57 25.51 37.61 42.12 7.97 ...

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)
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 効果量算出

# パッケージのインストール（初回のみ）
#install.packages("compute.es", dependencies = TRUE)
# パッケージの読み込み
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

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

d:「Cohenのd」
g:「Hedgesのg」

2.2 対応のないt検定の効果量

#mes(mean(a), mean(b), sd(a), sd(b), n.1, n.2)
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

2.3 対応なしの場合の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

2.4 対応のある場合の効果量d　

対応のない場合の数値を相関係数で補正する（p98参照）

# 事前テストと事後テストの相関係数は0.6487101
d.val <- res[, "d"] / sqrt(2 * (1 - 0.6487101))
# 95%信頼区間の計算
library("psych")
d.ci(d.val, n1 = 30)

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

3 ノンパラメトリック検定

3.1 Mann–Whitney

3.1.1 ファイルの読み込み

dat <- read.csv("data/data_ch6-1.csv", header = TRUE)

# パッケージのインストール（初回のみ）
#install.packages("coin", dependencies = TRUE)
# パッケージの読み込み
library("coin")

## Loading required package: survival

3.1.2 Mann–WhitneyのU検定

res <- wilcox_test(dat$score ~ factor(dat$class), distribution = "exact")
res

## 
##  Exact Wilcoxon-Mann-Whitney Test
## 
## data:  dat$score by factor(dat$class) (1, 2)
## Z = 2.3981, p-value = 0.01506
## alternative hypothesis: true mu is not equal to 0

# res@statistic@teststatisticは，上の結果のZ値と同じ
r <- abs(res@statistic@teststatistic) / sqrt(length(dat$score))
r

## [1] 0.489508

3.1.3 効果量rの95%信頼区間の算出

library("psych")
r.con(r, length(dat$score), p = .95, twotailed = TRUE)

## [1] 0.1072993 0.7456619

3.2 Wilcoxon

3.2.1 ファイルの読み込み

dat.2 <- read.csv("data/data_ch6-2.csv", header = TRUE)

# パッケージのインストール（初回のみ）
#install.packages("coin", dependencies = TRUE)
# パッケージの読み込み
library("exactRankTests")

##  Package 'exactRankTests' is no longer under development.
##  Please consider using package 'coin' instead.

## 
## Attaching package: 'exactRankTests'

## The following objects are masked from 'package:coin':
## 
##     dperm, pperm, qperm, rperm

3.2.2 Wilcoxonの符号付順位和検定

res.2 <- wilcox.exact(dat.2$pre, dat.2$post, paired = TRUE)
res.2

## 
##  Exact Wilcoxon signed rank test
## 
## data:  dat.2$pre and dat.2$post
## V = 74.5, p-value = 0.04126
## alternative hypothesis: true mu is not equal to 0

# 効果量rの計算
z <- qnorm(1 - (res.2$p.value / 2))
r <- z / sqrt(length(dat.2$pre) * 2)
r

## [1] 0.3726175

# 効果量rの95%信頼区間の算出
r.con(r, length(dat.2$pre * 2), p =.95, twotailed = TRUE)

## [1] -0.1725888  0.7430486

3.3 Kruskal-Wallis

3.3.1 ファイルの読み込み

dat.3 <- read.csv("data/data_ch6-3.csv", header = TRUE)

# Kruskal-Wallis検定
res.3 <- kruskal.test(dat.3$score ~ factor(dat.3$class))
res.3

## 
##  Kruskal-Wallis rank sum test
## 
## data:  dat.3$score by factor(dat.3$class)
## Kruskal-Wallis chi-squared = 10.879, df = 2, p-value = 0.004341

# 効果量rの計算
z.2 <- qnorm(1 - (res.3$p.value / 2))
r.2 <- abs(z.2) / sqrt(nrow(dat.3))
r.2

## [1] 0.4251917

# 効果量rの95%信頼区間の算出
r.con(r.2, nrow(dat.3), p =.95, twotailed = TRUE)

## [1] 0.1504330 0.6389761

3.4 多重比較

# 引数p.adjで"bonferroni"を指定すれば，Bonferroniの方法
pairwise.wilcox.test(dat.3[, 2], dat.3[, 1], p.adj = "holm", exact = FALSE, correct = FALSE)

## 
##  Pairwise comparisons using Wilcoxon rank sum test 
## 
## data:  dat.3[, 2] and dat.3[, 1] 
## 
##   1      2     
## 2 0.0022 -     
## 3 0.6329 0.0618
## 
## P value adjustment method: holm

3.5 Freedman

3.5.1 ファイルの読み込み

dat.4 <- read.csv("data/data_ch6-4.csv", header = TRUE)

# Freedman検定（データを行列に変換）
res.4 <- friedman.test(as.matrix(dat.4))
res.4

## 
##  Friedman rank sum test
## 
## data:  as.matrix(dat.4)
## Friedman chi-squared = 12.766, df = 2, p-value = 0.00169

# 効果量rの計算
z.3 <- qnorm(1 - (res.4$p.value / 2))
r.3 <- abs(z.3) / sqrt(nrow(dat.4))
r.3

## [1] 0.8708503

# 効果量rの95%信頼区間の算出
r.con(r.3, nrow(dat.4), p =.95, twotailed = TRUE)

## [1] 0.6149190 0.9608129

# 型の変更
dat.5 <- stack(dat.4) 
x <- dat.5[, 1]
y <- dat.5[, 2]
# 多重比較（引数pairedでTRUEを指定）
pairwise.wilcox.test(x, y, p.adj = "holm", exact = FALSE, paired = TRUE, correct = FALSE)

## 
##  Pairwise comparisons using Wilcoxon signed rank test 
## 
## data:  x and y 
## 
##         pre  post
## post    0.04 -   
## delayed 0.96 0.04
## 
## P value adjustment method: holm

6章_効果量

1 効果量

1.1 データ作成関数

1.1.1 データセットAの作成

1.1.2 データセットBの作成

1.1.3 t検定

1.2 因子型に変換

2 効果量算出

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

2.2 対応のないt検定の効果量

2.3 対応なしの場合のdの計算

2.4 対応のある場合の効果量d

3 ノンパラメトリック検定

3.1 Mann–Whitney

3.1.1 ファイルの読み込み

3.1.2 Mann–WhitneyのU検定

3.1.3 効果量rの95%信頼区間の算出

3.2 Wilcoxon

3.2.1 ファイルの読み込み

3.2.2 Wilcoxonの符号付順位和検定

3.3 Kruskal-Wallis

3.3.1 ファイルの読み込み

3.4 多重比較

3.5 Freedman

3.5.1 ファイルの読み込み

2.4 対応のある場合の効果量d　