裏 RjpWiki

文字通り,RjpWiki の裏を行きます
R プログラム コンピュータ・サイエンス 統計学

信頼区間が重なっていても,有意差があることもあるんだよ

2016年07月28日 | 統計学

「有意差なしのロゴです。有意差がないときにご利用ください」だって...

https://pbs.twimg.com/media/CoTKh__VYAEhjkT.jpg

図が傾いているし,左右の棒の基線がなんなのか?

一番問題なのは,エラーバーがなんなのか。

たとえば,左側の棒の下端が第一群の標本平均,右側の棒の上端が第二群の標本平均として,エラーバーはそれぞれの母平均の95%信頼区間を表しているものとする。

以下のテストデータで,提示されロゴが表す状態が再現できる。

n1 = 40; mean1 = -39; se1 =   6/qt(0.975, n1-1); sd1 = sqrt(n1)*se1
n2 = 30; mean2 =  83; se2 = 128/qt(0.975, n2-1); sd2 = sqrt(n2)*se2

x = scale(rnorm(n1))*sd1+mean1
y = scale(rnorm(n2))*sd2+mean2
d = data.frame(g=rep(1:2, c(n1, n2)), z=c(x, y))

plot(c(0, 3), c(-50, 220), type="n", xaxt="n", yaxt="n", xlab="", ylab="", bty="n")
w = 0.15
rect(1-w, c(0, mean1), 1+w, c(mean1, 0))
rect(2-w, c(0, mean2), 2+w, c(mean2, 0))
abline(h=0, col="red", lty=3)
arrows(1:2, c(mean1-6, mean2-128), 1:2, c(mean1+6, mean2+128), angle=90, length=w*0.7, code=3)
text(1+1.2*w, mean1, paste("mean =", mean1), pos=4)
text(2+1.2*w, mean2, paste("mean =", mean2), pos=4)
text(2+1.2*w, 0, "base line", pos=4)
text(1:2+1.2*w, mean1-se1*qt(0.975, n1-1), "mean-se = -45", pos=4)

では,このデータで平均値の差の検定を行うとどうなるか

> t.test(x, y, var.equal=TRUE)

    Two Sample t-test

data:  x and y
t = -2.2519, df = 68, p-value = 0.02756

となり,信頼区間が重なっているにもかかわらず「有意な差である」という結論になる。

もっとも,var.equal=FALSE にすると,サンプルサイズをいろいろ変化させても有意な差ではないという結果にはなる。

 

> t.test(x, y, var.equal=FALSE)

    Welch Two Sample t-test

data:  x and y
t = -1.9472, df = 29.13, p-value = 0.06121

いずれにせよ,「エラーバーが重なっているから有意差はない」とはいい切れないのである。

コメント

ダミー変数をカテゴリー変数に変換

2016年07月03日 | ブログラミング

奥村先生が,名古屋市のHPVVのデータの解析をしている

kaito = read.csv("kaito.csv", header=FALSE, colClasses="character", fileEncoding="UTF-8")
dim(kaito)
birth = ifelse(rowSums(kaito[,5:11]=="1") != 1, NA,
  ifelse(kaito[,5]=="1", 6,
  ifelse(kaito[,6]=="1", 7,
  ifelse(kaito[,7]=="1", 8,
  ifelse(kaito[,8]=="1", 9,
  ifelse(kaito[,9]=="1", 10,
  ifelse(kaito[,10]=="1", 11,
  ifelse(kaito[,11]=="1", 12, NA))))))))
table(birth, useNA="ifany")

birth はご本人も「もっとエレガントな方法があるだろうが,面倒なので…」と書いてあるが,簡単なので以下のように

birth2 = as.matrix(kaito[,5:11] == 1) %*% 6:12
birth2 = ifelse(birth2 < 6 | birth2 > 12, NA, birth2)
table(birth2, useNA="ifany")

kaito[,5:11] == 1 は,元のデータを colClasses="character" で読んだので,必要悪
掛けるベクトルの選択には注意(6:12 なら問題なし)
数値で読んでいればそのまま行列掛算
無視すべき(NA にすべき)データは,ifelse(birth2 < 6 | birth2 > 12, NA, birth2) で処理

それにしても,ダミー変数展開をそのままデータ入力するとは非効率甚だしい

---------

質問:以下の項目で当てはまるものにチェックしてください

□ Aです □ Bです □ Cです □ Dです ✓ Eです □ Fです □ Gです

0, 0, 0, 1, 0, 0, 0 って入力したのか...バカ

普通は

質問:以下の項目で当てはまる番号に○をつけてください

1. Aです 2. Bです 3. Cです 4. Dです 5. Eです 6. Fです 7. Gです 8. いずれでもない 9. わからない

一桁の数字で入力するだろ

また,いずれでもない,わからない などの項目も設けていないと,無回答なのか,該当しないのか,単に回答し忘れたのかなどすら区別できない

統計学以前,統計調査票(アンケート調査票)の作り方を知らない素人がやった調査なのか?

コメント