恣意的なカテゴリー分けで有意差を捏造する（黒い統計家シリーズ）
恣意的なカテゴリー分け
連続変数で有意にならなかったからカテゴリ変数にして有意にしてる論文を読んだという話を聞いた
これは面白い。つか、よくアクセプトされたな。
黒い統計家になった前提で、どの程度、この技が使えるのかやってみよう。
### 黒い統計家
nsim <- 1000
res <- matrix(0, nrow = nsim, ncol = 3)
for (j in 1:nsim) {
n <- 100
y <- rnorm(n, 10, 5)
age <- rnorm(n, 50, 15)
p.glm <- summary(glm(y ~ 1 + age, family = gaussian))$coefficients[2, 4]
y2 <- y[order(age)]
p.kuro <- numeric(n - 3)
for (i in 2:(n - 2)) {
p.kuro[i - 1] <- t.test(y2[1:i], y2[(i + 1):n], var.equal = FALSE)$p.value
}
res[j, ] <- c(p.glm, sum(p.kuro < 0.05), sum(p.kuro < p.glm))
}
res2 <- res[res[, 1] > 0.05, ]
summary(res2) #glmのP値 何個P<0.05のカテゴリー分けが可能か glmのP値より小さくなる回数
## V1 V2 V3
## Min. :0.0501 Min. : 0.00 Min. : 1
## 1st Qu.:0.2981 1st Qu.: 0.00 1st Qu.:31
## Median :0.5407 Median : 1.00 Median :47
## Mean :0.5324 Mean : 3.62 Mean :51
## 3rd Qu.:0.7760 3rd Qu.: 5.00 3rd Qu.:72
## Max. :0.9991 Max. :42.00 Max. :97
sum(res2[, 2] != 0)/nrow(res2) #glmでP>0.05になった場合に、カテゴリー分けを試行錯誤して、１つ以上のP<0.05の組合せを見つられる割合
## [1] 0.6033
なんと、普通の線形回帰で有意差が出なくっても、恣意的にカテゴリー分けを いろいろ試行錯誤すれば、約60％で有意差を捏造できる。 今回のシミュレーションでは、全く差がないデータでの話だから、恐ろしい までの確率である。
統計クイズ１（あなたの直感は正しいだろうか） t検定編
統計クイズ（あなたの直感は正しいだろうか）
＃ｔ検定の問題
各群10人ずつのデータが有る。
a群のアウトカムは平均10.7
b群のアウトカムは平均15.4であった。
アウトカムは数が大きいほど良いとする。
あなたはｂ群の方が良いんじゃないかと思っている。
これをｔ検定した所
## QUIZ
set.seed(1)
a <- rnorm(10, 10, 5.5)
b <- rnorm(10, 14, 5.5)
mean(a)
## [1] 10.73
mean(b)
## [1] 15.37
t.test(a, b, var.equal = FALSE)
##
## Welch Two Sample t-test
##
## data: a and b
## t = -2.015, df = 16.47, p-value = 0.06047
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -9.5122 0.2291
## sample estimates:
## mean of x mean of y
## 10.73 15.37
P=0.06であり、惜しくも有意差はなかった。
しかし、あと一息で有意差が出そうである。
あなたは黒い統計学者で、不適格のため脱落したがｂ群で 非常に良いアウトカム100を叩きだしたケースのことを思い出した。
これを,黒いあなたは事後的に組み入れようとした。
さて、新たなP値はいくらになるだろうか、Rを使わずに印象だけで 答えよ。
a群の平均値10.7
b群の平均値15.4
新たなデータを加えた新規b群の平均値23.1
選択肢は５つ
A:P=0.06
B:P=0.04
C:P=0.15
D:P=0.01
E:P<0.001
＃＃答えは下の方にスクロール＃＃＃
bb <- c(b, 100)
mean(bb)
## [1] 23.06
mean(bb) - mean(a)
## [1] 12.34
mean(b) - mean(a)
## [1] 4.642
t.test(a, bb, var.equal = FALSE)
##
## Welch Two Sample t-test
##
## data: a and bb
## t = -1.544, df = 10.59, p-value = 0.152
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -30.008 5.337
## sample estimates:
## mean of x mean of y
## 10.73 23.06
答えはCのP=0.15である。
平均値の差は拡大したのにP値はあがった。
100というのは超絶に外れ値である。 t検定は外れ値に対し脆弱である。分散の推定が崩壊するからだ。 超絶な外れ値の存在下では、どうやっても有意差が出ない。 結構Nが大きくても有意差は吹っ飛ぶので面白い。
はたして、印象だけで正しい値が想像できただろうか。
set.seed(1)
a <- rnorm(100, 10, 5.5)
b <- rnorm(100, 14, 5.5)
bb <- c(b, 1000)
mean(a)
## [1] 10.6
mean(b)
## [1] 13.79
mean(bb)
## [1] 23.56
t.test(a, b, var.equal = FALSE)
##
## Welch Two Sample t-test
##
## data: a and b
## t = -4.421, df = 197.2, p-value = 1.619e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -4.617 -1.769
## sample estimates:
## mean of x mean of y
## 10.60 13.79
t.test(a, bb, var.equal = FALSE)
##
## Welch Two Sample t-test
##
## data: a and bb
## t = -1.323, df = 100.5, p-value = 0.1887
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -32.381 6.466
## sample estimates:
## mean of x mean of y
## 10.60 23.56
上記の例でも、外れ値を入れると、平均値の差はむしろ拡大したのに 有意差がなくなってしまっている。
本日のトリビア t検定は外れ値に弱い！
