2015年3月17日火曜日

Rのtapply関数を使って、複数の度数分布表・ヒストグラムをまとめて出力する

tapplyを使えば、データフレームのカテゴリの列でグループ分けしつつ、値の列に対して関数を適用するということができます。

・・・と言っても分かりにくいと思いますので、具体的なサンプルをあげてみますと、

> # サンプルデータを作る
> 名前 <- c("A", "B",  "C", "D",  "E", "F")
> 性別 <- c("男", "男", "女", "男", "女", "女" )
> 身長 <- c(175,  165, 165,  170, 160, 155  )
> df <- data.frame(名前, 性別, 身長)
> df

  名前 性別 身長
1    A   男  175
2    B   男  165
3    C   女  165
4    D   男  170
5    E   女  160
6    F   女  155

>
> # ここから本題
> tapply(df$身長, df$性別, mean)
 
 女  男
160 170


性別でグループ分けして、身長の平均を算出する(mean関数を適用)という感じです。

第3引数にhist関数を指定してやれば、グループ分けを適用したあとに、複数の度数分布表を算出したり、複数の度数分布図(ヒストグラム)を書いたりすることもできます。

もう少し大きいサンプルデータを使いましょう。↓この本に載っていたサンプルを使わせていただきます。

データマイニング入門

↓こちらからダウンロードできる、"年収.csv"を使います。

データマイニング入門-Rで学ぶ最新データ解析 - 東京図書


↓こんな感じのサンプルデータです。

> d <- read.csv("年収.csv")
> head(d)
  学位 地位 性別 現職年数 学位年数  年収
1 博士 教授 男性       25       35 36350
2 博士 教授 男性       13       22 35350
3 博士 教授 男性       10       23 28200
4 博士 教授 女性        7       27 26775
5 修士 教授 男性       19       30 33696
6 博士 教授 男性       16       21 28516

tapplyでhistを適用しても、バババっとヒストグラムが次々描画される(そして消えていく)ことになってしまうので、ちょっと工夫が必要です。

工夫と言っても、戻り値を受け取って、所望の形式で出力させるだけですが。↓こんな感じ。

> # 度数の集計のため(描画が目的ではない)
> h <- tapply(d$年収, d$性別, hist, breaks=seq(15000, 40000, 5000))
>
> # 度数分布表の出力
> for ( i in 1:length(h) ) {
+   print ( names(h[i]) )
+   print ( data.frame(階級値=h[[i]]$mids, 度数=h[[i]]$counts) )
+   cat ( "\n") # 改行
+ }
[1] "女性"
  階級値 度数
1  17500    6
2  22500    5
3  27500    2
4  32500    0
5  37500    1

[1] "男性"
  階級値 度数
1  17500    9
2  22500   13
3  27500    9
4  32500    5
5  37500    2

それぞれに適用されたhistの戻り値がリストの形で h に格納されてますので、h[[1]], h[[2]], h[[3]], ... という感じで、それぞれにアクセスします。$midsには$breaksの中間の値が入っているので階級値として使えます。$countsには度数が入っています。

違うグループ分け(地位)でやってみましょう↓

> # 今度は、地位でグループ分け
> h <- tapply(d$年収, d$地位, hist, breaks=seq(15000, 40000, 5000))
>
> # 度数分布表の出力
> for ( i in 1:length(h) ) {
+   print ( names(h[i]) )
+   print ( data.frame(階級値=h[[i]]$mids, 度数=h[[i]]$counts) )
+   cat ( "\n") # 改行
+ }
[1] "教授"
  階級値 度数
1  17500    0
2  22500    2
3  27500   10
4  32500    5
5  37500    3

[1] "講師"
  階級値 度数
1  17500   15
2  22500    3
3  27500    0
4  32500    0
5  37500    0

[1] "助教"
  階級値 度数
1  17500    0
2  22500   13
3  27500    1
4  32500    0
5  37500    0

次は、ヒストグラムを書いてみましょう。前述のスクリプトでも描画自体はされていたのですが、複数なので上書きされてしまうのと、タイトルが入っていないのでどのグループのヒストグラムか分からない、というところをちょっと工夫します↓

# あらかじめグループ名を取得するために
h <- tapply(d$年収, d$性別, hist, breaks=seq(15000, 40000, 5000))
 
# カテゴリの数だけヒストグラムを縦に並べるための設定
par(mfrow=c(length(h),1))
 
# ヒストグラムの描画
for ( i in 1:length(h) ) {
  plot(h[[i]], main=names(h[i]), xlab="年収", col="gray")
}

tapplyでヒストグラムを描く(性別でグループ分け)

別のグループ分けを適用したいときは、tapplyの引数を1つ変えるだけ↓

# 今度は、地位でグループ分け
h <- tapply(d$年収, d$地位, hist, breaks=seq(15000, 40000, 5000))
 
# カテゴリの数だけヒストグラムを縦に並べるための設定
par(mfrow=c(length(h),1))
 
# ヒストグラムの描画
for ( i in 1:length(h) ) {
  plot(h[[i]], main=names(h[i]), xlab="年収", col="gray")
}

tapplyでヒストグラムを描く(地位でグループ分け)

カテゴリの種類が多い場合や、カテゴリが含む水準が多い場合は、いちいち手でやっていられないですからね。ちょっとExcelにはできない芸当だと思います。




2015年3月9日月曜日

片側検定/両側検定の選択が恣意的な気がする件をRのt検定で確認

こんな例を考えます。

10人の元々のテストの平均点は50点でした。点数をあげる施策を行った後、再度テストをしてみると、

after <- c(44, 47, 48, 51, 54, 56, 59, 62, 63, 65)

という点数が得られました。施策は効果があったかどうかを、有意水準5%で検定しましょう。

【Aさんの検定】

Aさんは施策に懐疑的でした。場合によっては施策の悪影響が出て、点数が下がる可能性もあると考えていますので、両側検定を行います。

> t.test(after, mu=50,  alternative="two.sided")

        One Sample t-test

data:  after
t = 2.1198, df = 9, p-value = 0.06306
alternative hypothesis: true mean is not equal to 50
95 percent confidence interval:
 49.67088 60.12912
sample estimates:
mean of x
     54.9


p値が0.05よりも大きいので、帰無仮説は棄却できず、「施策は効果があるとは言えない」と結論付けました。


【Bさんの検定】

Bさんは施策の効果を確信していました。点数が下がることはありえないと考えていますので、片側検定を行います。

> t.test(after, mu=50,  alternative="greater")

        One Sample t-test

data:  after
t = 2.1198, df = 9, p-value = 0.03153
alternative hypothesis: true mean is greater than 50
95 percent confidence interval:
 50.66264      Inf
sample estimates:
mean of x
     54.9


p値が0.05よりも小さいので、帰無仮説は棄却でき、「施策は効果があると言える」と結論付けました。


う~ん、気持ち悪い。

何が気持ち悪いって、効果に懐疑的だと厳しめの検定になって、効果に確信があると甘めの検定になるという関係です。

Aさんが施策反対派、Bさんが施策推進派だとすると、両者の対立の溝は埋められそうにないですね。

学術的にはこのへんの曖昧さの課題とかってないんでしょうか。統計の専門家に聞いてみたい。

だって、良かれと思ってやったことが悪影響を及ぼす、なんてことはよくあることでしょ。だとすると、片側検定していいケースってあんまりないような気がするんですが、どうなんでしょうか。




2015年3月6日金曜日

Rで対応のあるデータを比較する

例えば↓こんなデータがあったとします。40人のクラスがあって、1度目のテスト実施後に、なんらかの施策をして、その後、2度目のテストを実施したと。で、施策の効果は点数に表れているのか? みたいなのを想像していただけると理解しやすいかと。行名の 1, 2, ... , 40 はその学生の出席番号みたいなイメージで。

> df
   first second
1     56     75
2     50     59
3     53     61
4     60     62
5     45     59
6     51     60
7     48     70
8     55     79
9     50     67
10    44     50
11    50     64
12    56     49
13    42     53
14    52     53
15    57     41
16    39     63
17    50     54
18    45     54
19    54     59
20    48     74
21    48     58
22    47     59
23    47     57
24    43     54
25    56     76
26    54     39
27    48     59
28    45     43
29    50     59
30    51     52
31    50     58
32    48     70
33    52     61
34    51     50
35    43     62
36    44     39
37    49     64
38    47     45
39    52     68
40    46     64


で、いろいろな可視化をしてみようと思います。

1つ目は直線の傾きで見る感じです。matplot関数をつかえば直線が右肩上がりか、右肩下がりかで効果を把握できます。加えて、色分けすることで分かりやすくしてみました。

# 変化に応じて線の色を決める 
dif <-  df$second - df$first
color <- ifelse( dif > 0, "blue",             # 上がっていれば青
           ifelse( dif < 0, "red", "black") ) # 下がっていれば赤、変化なしは黒
matplot(t(df), type="l", col=color, lty=1, xaxt="n", ylab="score")
axis(1, c(1,2), c("first", "second"))
matplot(df, xlab="学生番号", ylab="点数")

matplotの例

全体的に上がっている人が多そうなことが分かります。また、ばらつきが大きくなっていることも見てとれますね。

t()関数で転置させているのは、first, secondを横軸方向に取るためです。転置させずにmatplotすると、横軸方向は学生番号(1~40)になります↓

matplot(df, xlab="学生番号", ylab="点数")
matplotの例

これはこれで、ありのような気がしますが、1人の学生に着目しての、上がった/下がったというのは分かりにくいですね。縦方向に照合しにくいとか、差分の距離が分かりにくいとか、そのへんが原因でしょうか。

差分に対してヒストグラムを書いてみるというのも、1つのアプローチですね。

hist( dif, breaks=seq(-50, 50, 10), col=c(rep("red",5), rep("blue",5)) )


差の値でヒストグラムを描く

上がった人が多いことや、どの程度上がった/下がったのかという分布が一目で分かりますね。

変り種としては、1回目を横軸に、2回目を縦軸に取って散布図を描いてみるというものです。

plot(df$first, df$second, asp=1, xlim=c(0,100), ylim=c(0,100))
abline( 20, 1, col="blue", lty=1 )
abline( 10, 1, col="blue", lty=2 )
abline(  0, 1 )
abline(-10, 1, col="red" , lty=2 )
abline(-20, 1, col="red" , lty=1 )
 
legend( "topleft", c("20点UPライン", "10点UPライン", "10点DOWNライン", "20点DOWNライン"),
        col=c("blue", "blue", "red", "red"), lty=c(1, 2, 2, 1))

前をx軸、後をy軸にとって散布図を描く

y = x の直線(黒の実線)の上に乗っていれば変化なし、上側なら点数UP、下側なら点数DOWNという感じです。直線からの離れ具合が、UP/DOWNの程度に対応しているイメージですね。

「誰が」というところを知りたいなら、該当する学生番号をtext関数でプロットする手もあります。

plot(df$first, df$second, asp=1, xlim=c(0,100), ylim=c(0,100), type="n")
text(df$first, df$second, row.names(df))

前をx軸、後をy軸にとって散布図を描く(text関数でIDをプロット)

ごちゃごちゃして見えにくい場合は、字を小さくするなり、画像のサイズを大きくするなりしてください。

出力先をPDFにし、ごく小さい字でプロットして、拡大・縮小しながら見るというのも結構いい感じです。PDFに出力する場合は↓下記が参考になると思います。

Rでshapeを出力するときはPDFファイルにすると便利 - Rプログラミングの小ネタ


shapeに限らず、こまごましたものをプロットするときはPDFがおすすめです。

最後のサンプルは、個々の学生の上がり/下がり、その幅に着目したいような場合の例です。色(colorという変数に格納)は、1つ目のサンプルで設定したものを使いまわしています。

plot(1:n, df$first, xlim=c(0,40), ylim=c(0,100), xlab="学生番号", ylab="点数", col=color)
arrows(1:n, df$first, 1:n, df$second, length=0.1, col=color)

増減をarrowsで表現
これだと、変化幅が良く分かりますね。

・・・などという感じで、いろいろ紹介してみました。




Rでエラーバー付きのプロットを行う

例えば↓こんなサンプルデータがあったとしましょう。属するグループと観測値の組みたいな感じですね。

> x1

   group value
1      A  2.66
2      A  2.24
3      A  2.41
4      A  1.12
5      A  2.48
6      A  2.33
7      A  2.90
8      A  4.35
9      A  3.82
10     A  3.45
11     B  6.65
12     B  4.61
13     B  5.26
14     B  5.24
15     B  7.26
16     B  5.97
17     B  2.90
18     B  2.37
19     B  4.94
20     B  4.08
21     C  5.63
22     C  5.75
23     C  4.87
24     C  8.14
25     C  1.01
26     C  1.58
27     C  3.76
28     C  8.12
29     C  5.95
30     C  4.48


で、グループごとの平均やばらつきにどんな傾向があるか見るために、エラーバー付きでプロットしたいような状況だと考えてください。

こういうときは、gplotsパッケージのplotmeans関数が便利です。

install.packages("gplots")
library(gplots)
plotmeans(value~group, data=x1, connect=F, ylim=c(0,10))
plotmeansによるエラーバー付きプロット

ところで、エラーバーの長さは何を表しているのでしょうか?

> ?plotmeans

でマニュアルの引数のところを見ると、

p     confidence level for error bars. Defaults to 0.95.

と書いてあります。p= で指定することによって、エラーバーの信頼水準を指定することができ、そのデフォルト値は0.95であると書いてあるので、つまりエラーバーは信頼区間を表しているということですね。

確認してみましょう。

グループAだけ抽出して、t.test関数を使ってt検定をし、その95パーセント信頼区間を求めてみます。

> tt <- t.test(x1[x1$group=="A",]$value)
> tt

        One Sample t-test

data:  x1[x1$group == "A", ]$value
t = 9.6159, df = 9, p-value = 4.952e-06
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 2.122944  3.429056
sample estimates:
mean of x
    2.776

2.122944 ~ 3.429056 という結果が出ました。これを先ほどのplotmeansの結果に追記してみます。

abline(h=tt$conf.int[1], lty=2)
abline(h=tt$conf.int[2], lty=2)
plotmeansのエラーバーは信頼区間であることを確認

うん、見事に一致しました。

ということで、plotmeansにおけるエラーバーはt検定における信頼区間を現していて、デフォルトは0.95だけど、p= で指定できるということが分かってもらえたと思います。

で、以降は補足です。

データの形が違うってことがありますよね。例えば↓こんな感じ。

> x2
      A    B    C
1  2.66 6.65 5.63
2  2.24 4.61 5.75
3  2.41 5.26 4.87
4  1.12 5.24 8.14
5  2.48 7.26 1.01
6  2.33 5.97 1.58
7  2.90 2.90 3.76
8  4.35 2.37 8.12
9  3.82 4.94 5.95
10 3.45 4.08 4.48

こんな場合は、plotmeansの入力に合わせて、形を組み替えちゃえばいいと思います。

> # plotmeansの入力データの形式に直す
> x3 <- data.frame(value=unlist(x2), group=rep(names(x2), each=nrow(x2)))
> x3
    value group
A1   2.66     A
A2   2.24     A
A3   2.41     A
A4   1.12     A
A5   2.48     A
A6   2.33     A
A7   2.90     A
A8   4.35     A
A9   3.82     A
A10  3.45     A
B1   6.65     B
B2   4.61     B
B3   5.26     B
B4   5.24     B
B5   7.26     B
B6   5.97     B
B7   2.90     B
B8   2.37     B
B9   4.94     B
B10  4.08     B
C1   5.63     C
C2   5.75     C
C3   4.87     C
C4   8.14     C
C5   1.01     C
C6   1.58     C
C7   3.76     C
C8   8.12     C
C9   5.95     C
C10  4.48     C

行名が気に入らない場合は、↓こんな感じで消してしまいましょう。

> row.names(x3) <- NULL
> x3
   value group
1   2.66     A
2   2.24     A
3   2.41     A
4   1.12     A
5   2.48     A
6   2.33     A
7   2.90     A
8   4.35     A
9   3.82     A
10  3.45     A
11  6.65     B
12  4.61     B
13  5.26     B
14  5.24     B
15  7.26     B
16  5.97     B
17  2.90     B
18  2.37     B
19  4.94     B
20  4.08     B
21  5.63     C
22  5.75     C
23  4.87     C
24  8.14     C
25  1.01     C
26  1.58     C
27  3.76     C
28  8.12     C
29  5.95     C
30  4.48     C

以下は、さらに蛇足。

Rでエラーバーを表示させる方法といえば、私は↓この方法しか知らなかったんですね。

誤差範囲 | Rで棒グラフと折れ線グラフにエラーバーを付ける | バイオスタティスティクス


矢印を書くarrows関数というのがあるんですが、矢尻の角度を矢の本体の棒に対して90度にすることで、エラーバーっぽく見せるというやり方です。矢の長さを決めるのに必要な値は自分で事前に計算しておくという感じで、かなり手作り感覚あふれる描画方法です。

で、Rにはエラーバーを表示する関数はないのかなあと思っていたところで、↓この書籍で、gplotsパッケージとplotmeans関数の存在を知りました。

Rによる統計的検定と推定
Rによる統計的検定と推定

ちなみに、この書籍の良いところは、それぞれの手法に対して、Rによる検定コマンドの紹介の前に、ちゃんと可視化の例が載っているところですね。

検定の前にちゃんとデータの可視化を行っている
「Rによる統計的検討と推定」より

上記は、2つの平均値(対応あり)の検定を紹介しているページですが、t.testによる検定の前に、matplotによる可視化でデータの傾向を確認しています。可視化をすっ飛ばして、いきなり検定するとかありえないですもんね。

話を戻すと、前述のリンク先(バイオスタティスティクス)のarrows関数による方法は、それでもやっぱり有用だったりします。

というのも、plotmeansによるエラーバーは信頼区間を現していましたが、エラーバーの長さって他の値で決めることもあるじゃないですか。1σとか。plotmeansではそれはできないのですが、arrowsを使う方法ならば、自分で値を計算するから、任意の指標でのエラーバー付けができますからね。これはこれでいろいろ使えて便利だと思います。