日誌

Rメモ(tidyverse) >> 記事詳細

2017/09/09

janitorパッケージでtidyな度数分布表やクロス表を作成

Tweet ThisSend to Facebook | by mtsuchiya
カテゴリ変数の水準ごとの人数,また変数同士のクロス表を算出する場面は多いと思います。そうした際に,tidyverseの影響を受けたtidyverse-orientedなパッケージ,janitorが便利だったので,使い方を整理してみました。

使用するデータは,psychパッケージのbfiデータです。データについて詳しくはこちら参照。

それではさっそくパッケージの読み込みからです。

library(tidyverse)
library(janitor)
library(psych)

中の変数確認。

names(bfi)
 
[1] "A1"        "A2"        "A3"        "A4"      
 [5] "A5"        "C1"        "C2"        "C3"      
 [9] "C4"        "C5"        "E1"        "E2"      
[13] "E3"        "E4"        "E5"        "N1"      
[17] "N2"        "N3"        "N4"        "N5"      
[21] "O1"        "O2"        "O3"        "O4"      
[25] "O5"        "gender"    "education" "age"


1.度数分布(欠損値がない場合)


度数を算出するには当該変数にtable()関数を適用するのがよく用いられていると思います。性別のgender変数で確認してみます。

table(bfi$gender)

   1    2
 919 1881

割合を算出するには,さらにprop.table()で囲む

prop.table(table(bfi$gender))


        1         2
0.3282143 0.6717857

dplyrのcount( )関数でも同じことができます。

bfi %>%  count(gender)

# A tibble: 2 x 2
  gender     n
   <int> <int>
1      1   919
2      2  1881

%を出したい場合は,mutateでそれぞれの水準のnを合計nで割った式を書き,新しくfreq変数を作成します。でもちょっとコードが長く感じますよね。

bfi %>%
  count(gender) %>%
  mutate(freq = n/sum(n))

# A tibble: 2 x 3
  gender     n      freq
   <int> <int>     <dbl>
1      1   919 0.3282143
2      2  1881 0.6717857

そこで,便利なのがjanitorパッケージのtabyl( )関数。これなら%>%でつないで関数に変数名を入れるだけでOK

bfi %>% tabyl(gender)

  gender    n   percent
1      1  919 0.3282143
2      2 1881 0.6717857


2.度数分布(欠損値がある場合)


では,次のような場合はどうでしょう。教育歴を表す変数,educationです。

table(bfi$education)

   1    2    3    4    5 
 224  292 1249  394  418

実は,この変数には欠損値があるので,それを表示させる引数を追加します。

table(bfi$education, useNA = "ifany")


   1    2    3    4    5 <NA>
 224  292 1249  394  418  223

割合はさらに囲む

prop.table(table(bfi$education, useNA = "ifany"))


         1          2          3          4          5
0.08000000 0.10428571 0.44607143 0.14071429 0.14928571
      <NA>
0.07964286


なんだかごちゃごちゃしてて見にくいですね。

dplyrのcountを使えば,NAはデフォルトで出してくれます。しかも結果もデータフレームで扱いやすいです。

bfi %>% count(education)

# A tibble: 6 x 2
  education     n
      <int> <int>
1         1   224
2         2   292
3         3  1249
4         4   394
5         5   418
6        NA   223

割合を出すには,mutateでnを合計で割った式を書きます

bfi %>%
   count(education) %>%
   mutate(freq = n/sum(n))

  education     n       freq
      <int> <int>      <dbl>
1         1   224 0.08000000
2         2   292 0.10428571
3         3  1249 0.44607143
4         4   394 0.14071429
5         5   418 0.14928571
6        NA   223 0.07964286

欠損値を無視したい場合はfilterで指定

bfi %>%
  filter(!is.na(education)) %>%
  count(education) %>%
  mutate(freq=n/sum(n))


  education     n       freq
      <int> <int>      <dbl>
1         1   224 0.08692278
2         2   292 0.11331005
3         3  1249 0.48467210
4         4   394 0.15289096
5         5   418 0.16220411


それではtabyl( )の便利さを見てみましょう

bfi %>% tabyl(education)


  education    n    percent valid_percent
1         1  224 0.08000000    0.08692278
2         2  292 0.10428571    0.11331005
3         3 1249 0.44607143    0.48467210
4         4  394 0.14071429    0.15289096
5         5  418 0.14928571    0.16220411
6        NA  223 0.07964286            NA


基本的なコードで欠損値を削除した時の値なども出してくれる。
さらに,総計の数が欲しいときはadorn_totals()関数を加えます。

bfi %>%
  tabyl(education) %>%
  adorn_totals(which = "row")


  education    n    percent valid_percent
1         1  224 0.08000000    0.08692278
2         2  292 0.10428571    0.11331005
3         3 1249 0.44607143    0.48467210
4         4  394 0.14071429    0.15289096
5         5  418 0.14928571    0.16220411
6      <NA>  223 0.07964286            NA
7     Total 2800 1.00000000    1.00000000



3.クロス表


次はgenderとeducationでクロス表をかいてみます。まずはtable( )で。


table(bfi$education,bfi$gender)
  
          1   2
      1  93 131
      2 103 189
      3 356 893
      4 134 260
      5 152 266  


割合を出してみます。多くの場合,興味あるのはrowの%かと思います。
例えば,教育歴が1の人達では女性が○○%,といった具合です。

prop.table(table(bfi$education,bfi$gender))

             1          2
  1 0.03608847 0.05083430
  2 0.03996896 0.07334109
  3 0.13814513 0.34652697
  4 0.05199845 0.10089251
  5 0.05898331 0.10322080


数字が出てきましたが,これは,全体から各セルの割合を見た数字なので,欲しい数値ではありません。janitorでは,crosstab( )という便利関数があります。

bfi %>% crosstab(education,gender)


  education   1   2
1         1  93 131
2         2 103 189
3         3 356 893
4         4 134 260
5         5 152 266
6        NA  81 142

rowの%を見たい時は,オプションにpercent = "row"を追加します。

bfi %>% crosstab(education,gender, percent = "row")

  education         1         2
1         1 0.4151786 0.5848214
2         2 0.3527397 0.6472603
3         3 0.2850280 0.7149720
4         4 0.3401015 0.6598985
5         5 0.3636364 0.6363636
6        NA 0.3632287 0.6367713


このように欲しい数値が出てきました。
prop.table()と同じく,全体に対する各セルの割合が欲しい場合は,引数にそう指定します。

bfi %>% crosstab(education,gender, percent = "all")


  education          1          2
1         1 0.03321429 0.04678571
2         2 0.03678571 0.06750000
3         3 0.12714286 0.31892857
4         4 0.04785714 0.09285714
5         5 0.05428571 0.09500000
6        NA 0.02892857 0.05071429


NAを除外したい場合は,引数 show_na = FALSE を加えます

bfi %>% crosstab(education,gender, percent = "all", show_na = FALSE)

  education          1          2
1         1 0.03608847 0.05083430
2         2 0.03996896 0.07334109
3         3 0.13814513 0.34652697
4         4 0.05199845 0.10089251
5         5 0.05898331 0.10322080

さらに,adorn_crosstab( )関数を使うとそのまま表にできる感じの出力もできます

bfi %>%
  crosstab(education,gender, show_na = FALSE) %>%
  adorn_crosstab()


  education           1           2
1         1 41.5%  (93) 58.5% (131)
2         2 35.3% (103) 64.7% (189)
3         3 28.5% (356) 71.5% (893)
4         4 34.0% (134) 66.0% (260)
5         5 36.4% (152) 63.6% (266)

22:39 | 投票する | 投票数(0) | コメント(0)