Rプログラミングの小ネタ

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 <- rea...

こんな例を考えます。 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    ...

例えば↓こんなデータがあったとします。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     ...

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

「シミュレーション」なんていうと、難しそうな気がしてしまいますが、今回のは「乱数を使って試してみる」くらいのものです。 例えば、10万人の学生がいます。彼らの身長は正規分布に従っているとします。また、彼らのテストの点数も正規分布に従っているとします。 このとき、「身長＋点数」という確率変数を考えると、これは正規分布に従っているでしょうか？ これを数学的に証明するのではなく、乱数を使ったシミュレーションで、「身長＋点数もきっと正規分布に従うんじゃないか」ということを確認しようというわけです。 まずは↓こんな感じでデータを生成しておきます。 # 身長 # 平均 ：170cm # 標準偏差： 5cm heights <- rnorm ( 100000 , 170 , 5 )   # 点数 # 平均 ： 50 # 標準偏差： 10 scores <- rnorm ( 100000 , 50 , 10 ) ↓身長と点数のそれぞれでヒストグラムを描いてみましょう。 brk <- seq ( 0 , 300 , 1 ) # ヒストグラムの刻み hist ( heights , breaks=brk , ylim= c ( 0 , 10000 ) ) 正規分布に従う「身長」のヒストグラム hist ( scores , breaks=brk , ylim= c ( 0 , 10000 ) ) 正規分布に従う「点数」のヒストグラム ↓両者を足してヒストグラムを描いてみると、 hist ( heights + scores , breaks=brk , ylim= c ( 0 , 10000 ) ) 「身長＋点数」は正規分布に従うか？ どうやら正規分布になっているっぽいです。 では、どんな正規分布（平均、分散）になっているのでしょうか。 > mean(heights) [1] 169.9721 > mean(scores) [1] 50.04487 > mean(heights + scores) [1] 220.017 「身長...

例えば、↓こんなヒストグラムがあったとしましょう。 普通にhist関数で描画したもの。 30～40、40～50の階級の度数はゼロ？ ぱっと見、30～40と40～50の階級の度数はゼロなんだな、と見えちゃいますが、そうではありません。 ↓実はこんなスクリプトでサンプルデータを作ったので、 d <- c ( rep ( 5 , 10000 ) , rep ( 15 , 1000 ) , rep ( 25 , 100 ) , rep ( 35 , 10 ) , rep ( 45 , 1 ) ) hist ( d , breaks= seq ( 0 , 50 , 10 ) ) 30～40の階級には10の度数、40～50の階級には1の度数があるんですよね。でも、0～10の階級の10000の度数が多すぎるために、つぶれて見えなくなってしまったという状態。 こんなときは、棒グラフの上のところに度数の数字を表示させてやれば、少ない度数の階級についても確認ができます。 hist ( d , breaks= seq ( 0 , 50 , 10 ) , labels =T ) hist関数にlabels=Tオプションを指定し、度数を表示させたもの labels=T オプションを指定するだけなので、とても簡単です。 が、実は私、このオプションの存在を知らず、いままで↓下記のような方法で書いてました。（このへんの紆余曲折が、この記事を書こうと思い立った理由だったりします） h <- hist ( d , breaks= seq ( 0 , 50 , 10 ) ) text ( h$mids , h$counts , h$counts ) hist関数の戻り値を、text関数を使ってプロット hist関数の戻り値（変数hに格納してます）には、階級値や度数の情報が入っているので、それをtext関数で表示させるというもの。 文字と棒グラフのborderが重なって見にくいと思う場合は、↓こうやればOK。 h <- hist ( d , breaks= seq ( 0 , 50 , 10 ) , ylim= c ( 0 , 10100 ) ) ...

duplicated関数の引数にデータフレームを指定すると、すべての項目（列）を組み合わせた上での重複チェックができますよ、という話を以前書きました↓ Rのデータフレームで複数項目を組み合わせた重複のチェック - Rプログラミングの小ネタ また、上記の記事では、すべてではなく一部の項目のみを組み合わせ対象にしたい場合は、interaction関数を使って新たにチェック用の列を作ってはどう？、というやり方も紹介しました。 が、duplicated関数がすべての列を対象にする動きなら、対象にしたい列だけで構成された新たなデータフレームを生成して、duplicated関数に渡した方がシンプルではないか、ということに気づきました。 例えば↓こんなデータがあって、 > d    名字 名前 年齢 1  伊東 二郎   52 2  佐藤 五郎   52 3  遠藤 太郎   53 4  後藤 太郎   53 5  近藤 五郎   28 6  遠藤 二郎   33 7  後藤 四郎   52 8  須藤 次郎   28 9  伊東 二郎   35 10 遠藤 太郎   23 11 遠藤 史郎   38 12 工藤 吾郎   53 13 須藤 太郎   22 14 武藤 二郎   28 15 工藤 吾郎   53 16 伊東 二郎   22 17 遠藤 四郎   33 18 後藤 史郎   38 19 加藤 一郎   28 20 加藤 二郎   22 名字と名前だけを対象として重複チェックをしたい場合は、 > sum(duplicated( data.frame(d$名字, d$名前) )) [1] 4 とか、 > d[duplicated( data.frame(...