【R】e-Statのデータを加工してggplotで作図


だんだんRを使い慣れてきたので、アウトプットがてらちょっとやってみたかったことをまとめてみました。
Qiita初投稿なので、手探りですがどうぞよろしくお願いします。

追記：タイトル変更しました
追記２：無事twitterのOGPも更新できました！（参考 https://qiita.com/AumyF/items/5944a27b3a41036040a2 ）thx!

# e-Statに公開されているデータを使ってなんかしたい
「e-Statに公開されているデータを使ってなんかしてみたい…」
そう思ったことのある人はきっとわたし意外にもいるはず…。
いろいろな情報が公開されているのは知ってるけど、前処理がいろいろと複雑そうだし、何をどう使えるのかが未知数…ってのがわたしの個人的な感想です。
習うより慣れろ、ってことで何かしらやってみたいと思います！

e-Statで公開されているデータの中で、特にわたしが興味を持っているのが「国民健康・栄養調査」のデータなんですよねぇ…。
例えば、日本人がたんぱく質とか脂質とかの栄養素を、どの食品からたくさん食べているのかって知りたくないですか？
知りたいよね！（強引）
せっかくなので、グラフになってたら見やすそうですよね。
よし、それつくる。

というわけで、今日は**日本人がたんぱく質をどの食品からたくさん食べているのか、Rを使って調べてみたいと思います！**
ついでに脂質とか炭水化物もみてみよかな。
棒グラフで横にずらっと食品の名前が並んでたらおもしろいのではなかろうか。

# 環境と前提条件
今回使うOSとRのバージョンは以下の通りです

- macOS Catalina 10.15.7
- R version 4.0.2 (2020-06-22)
- R Studio Version 1.3.1093

今回はとりあえず知ってる知識でなんかやってみたい！！を叶えようと思います。
パッケージは`tidyverse`と`magrittr`を使う予定です。パイプ`%>%`とかtidyverseの中に含まれている各パッケージについての説明は割愛しますのでご了承ください。
後々、これは説明が必要欲しいなぁってところが出てきたら、随時追記して行く予定です。

APIを使う方法はまだ勉強中なので、ひとまずcsvデータをDLする方法でやってみます。


# 1.e-Statから国民健康・栄養調査のデータをDLする
とりあえず下記のURLからe-Statにアクセス。
https://www.e-stat.go.jp/

キーワード検索から「国民健康・栄養調査」を検索します。
栄養素ごとに、どの食品からどれくらいの量を食べているのかがわかるデータ（食品群別栄養素等摂取量）を手に入れようと思います。
検索結果はひとつしか出てこないので、表示されるままクリックしてここまできます。


![スクリーンショット 2021-01-06 16.22.11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1005648/e4350b72-2223-9456-5bf5-f4d710f3b03b.png)

**データベース**をクリックして　**年次**　→　**2017年**　→　**表番号９のDB**（画像赤枠） の順に進みます。
![スクリーンショット 2021-01-06 17.12.39.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1005648/a08defb5-2844-ef1f-745d-ef486f7d2d04.png)



この画面が表示されたら、特になにもいじらずそのまま右の**ダウンロード**をクリック。
![スクリーンショット 2021-01-06 16.24.16.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1005648/ce748422-1a4f-5071-0d96-345ab0313ede.png)

表ダウンロードのページに飛ぶので、ヘッダの出力：<b>出力しない</b>、コードの出力：<b>出力しない</b>　を選んで（他はそのまま）、下部に表示されるダウンロードボタンを押します。
お好みの場所、ファイル名で保存します。
![スクリーンショット 2021-01-06 17.31.15.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1005648/695a9f00-bfb5-c6bb-7247-0e9f13993dc5.png)

データがDLできたら、ファイルを開いて中身を確認してみます。
>Macだと勝手にNumbersで開いちゃうけど、ここはExcelで開きます

![スクリーンショット 2021-01-06 17.32.17.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1005648/95768491-7cb4-04c7-bfc7-680e8c5721e4.png)
ちゃんとA1からデータで埋まっています。やったぜ
これで必要なデータが入手できました！



# 2.DLしたデータの読み込み

DLしたデータをRに読み込んでいきましょう。まずは`tidyverse`と`magrittr`を召喚。

```r:召喚の呪文
library(tidyverse)
library(magrittr)
```

わたしはlibraryを2回呼び出すのが面倒なので、`pacman`というパッケージで一気に召喚してます。

```r:さくっとpacmanで召喚するパターン
# pacmanが入ってなかったらインストール
if(!(require(pacman))) install.packages("pacman")

# tidyverseとmagrittrを呼び出す
pacman::p_load(pacman, tidyverse, magrittr)
```

`readr()`でDLしたデータを読み込むと…

```r:データの読み込み
data <- read_r('ファイル名.csv')

> data <- read_csv('ファイル名.csv')
 make.names(x) でエラー: 
   '<94>N<97><ee><8a><4b>��' に不正なマルチバイト文字があります 
```

ってな感じでエラーが出ます。
不正なマルチバイト文字ってことなので、文字のエンコードを指定する必要がありそうです。
日本語のデータなので、Shift-JISに指定してみます。

```r:Shift-JISでデータを読み込む
data <- read_csv('ファイル名', locale = locale(encoding = "SHIFT-JIS"))

> data <- read_csv('INPUT/FEH_00450171_210106173118.csv', locale = locale(encoding = "SHIFT-JIS"))
Parsed with column specification:
cols(
  .default = col_double(),
  年齢階級 = col_character(),
  `時間軸(年次)` = col_character(),
  食品群 = col_character(),
  `/栄養素等` = col_logical(),
  `摂取量【g】` = col_number(),
  `エネルギー【kcal】` = col_number(),
  `ナトリウム【mg】` = col_number(),
  `カリウム【mg】` = col_number()
)
See spec(...) for full column specifications.
```

Environmentにdataが追加されたら読み込みOK！
各列の変数型は`readr()`が適宜決めてくれたみたいです。これじゃねぇぞ、ってときは指定することもできますが、今回は問題なさそうなのでそのままで。

> ちなみにこれ、メニューバーから File > Import Dataset > From Text(readr)... ってやっても同じことができます。
> プレビュー画面が見れるので、探り探り読み込む時は、最初にそこから確認することが多いかも。


# 補足：ヘッダとコードを出力した場合の読み込み
もし、デフォルトのままヘッダとコードを出力していると以下のようにcsv出力されているはずです。
![スクリーンショット 2021-01-06 16.39.56.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1005648/2972ff98-f9d2-e8b0-90ef-87cbd3e2850d.png)
ヘッダーを出力していると、読み込みの際にエラーが出るので1行目から11行目をスキップする必要があります。
`skip = `を追加してどこまでスキップする行数を指定すればOKです。

```r:データの読み込み（ヘッダが含まれている場合）
# 11行目まではスキップして12行目から読み込んでね
data <- read_csv('ファイル名.csv', locale = locale(encoding = "SHIFT-JIS"), skip = 11) 
```

# 3.読み込んだデータをざくっと眺める

読み込んだデータの中身を`str()`で確認してみます。

```r:中身の確認
data %>% str()

> data %>% str()
tibble [262 × 38] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ 年齢階級                    : chr [1:262] "20歳以上" "20歳以上" "20歳以上" "20歳以上" ...
 $ 時間軸(年次)                : chr [1:262] "2017年" "2017年" "2017年" "2017年" ...
 $ 食品群                      : chr [1:262] "総量" "動物性食品" "植物性食品" "穀類" ...
 $ /栄養素等                   : logi [1:262] NA NA NA NA NA NA ...
 $ 摂取量【g】                 : num [1:262] 2108 320 1788 424 305 ...
 $ エネルギー【kcal】          : num [1:262] 1914 471 1443 760 513 ...
 $ たんぱく質【g】             : num [1:262] 70.5 37.8 32.7 15.2 7.7 7.5 0.1 7 0.3 3.3 ...
 $ 脂質【g】                   : num [1:262] 58.9 29 29.9 4.7 1 0.9 0 3.6 0.1 1.8 ...
 
# 長かったので以降は割愛
```

年齢階級と時間軸（年次）って列が気になりました。`summary.factor()`を使ってさらに確認してみます。

```r:年齢階級の要約
data$年齢階級 %>% summary.factor()

> data$年齢階級 %>% summary.factor()
20歳以上 
     262 
```
年齢階級の列は同じ情報だけのようなので特に気にする必要はなさそうです。

```r:時間軸（年次）の要約
data$`時間軸(年次)` %>% summary.factor()

> data$`時間軸(年次)` %>% summary.factor()
2016年 2017年 
   131    131 
```
あれ？2017年のデータをDLしたつもりが、2016年の情報も入ってるみたいですね。
今回はとりあえず、2017年のデータだけを使いましょう。

`slice()`を使って最初の数行の内容も確認してみます。

```r:最初の数行だけ覗いてみる
data %>% slice()

> data %>% slice()
# A tibble: 262 x 38
   年齢階級 `時間軸(年次)` 食品群 `/栄養素等` `摂取量【g】` `エネルギー【kcal】`… `たんぱく質【g】`… `脂質【g】`
   <chr>    <chr>          <chr>  <lgl>               <dbl>            <dbl>            <dbl>       <dbl>
 1 20歳以上 2017年         総量   NA                 2108.              1914             70.5        58.9
 2 20歳以上 2017年         動物性食品… NA                  320.               471             37.8        29  
 3 20歳以上 2017年         植物性食品… NA                 1788               1443             32.7        29.9
 4 20歳以上 2017年         穀類   NA                  424.               760             15.2         4.7
 5 20歳以上 2017年         米・加工品… NA                  305.               513              7.7         1  
 6 20歳以上 2017年         米     NA                  302                506              7.5         0.9
 7 20歳以上 2017年         米加工品… NA                    3.4                7              0.1         0  
 8 20歳以上 2017年         小麦・加工… NA                  108.               229              7           3.6
 9 20歳以上 2017年         小麦粉類… NA                    3.7               14              0.3         0.1
10 20歳以上 2017年         パン類（菓… NA                   36                 98              3.3         1.8
# … with 252 more rows, and 30 more variables: `脂肪酸_飽和脂肪酸【g】` <dbl>,
#   `脂肪酸_一価不飽和脂肪酸【g】` <dbl>, `脂肪酸_n-6系脂肪酸【g】` <dbl>, `脂肪酸_n-3系脂肪酸【g】` <dbl>,
#   `コレステロール【mg】` <dbl>, `炭水化物【g】` <dbl>, `食物繊維_総量【g】` <dbl>,
#   `食物繊維_水溶性【g】` <dbl>, `食物繊維_不溶性【g】` <dbl>, `ビタミンＡ_RE【μｇＲＥ】` <dbl>,
#   `ビタミンＤ【μｇ】` <dbl>, `ビタミンＥ【mg】` <dbl>, `ビタミンＫ【μｇ】` <dbl>,
#   `ビタミンＢ1【mg】` <dbl>, `ビタミンＢ2【mg】` <dbl>, `ナイアシン【mgNE】` <dbl>,
#   `ビタミンＢ6【mg】` <dbl>, `ビタミンＢ12【μｇ】` <dbl>, `葉酸【μｇ】` <dbl>,
#   `パントテン酸【mg】` <dbl>, `ビタミンＣ【mg】` <dbl>, `ナトリウム【mg】` <dbl>,
#   `食塩相当量_g【g】` <dbl>, `カリウム【mg】` <dbl>, `カルシウム【mg】` <dbl>, `マグネシウム【mg】` <dbl>,
#   `リン【mg】` <dbl>, `鉄【mg】` <dbl>, `亜鉛【mg】` <dbl>, `銅【mg】` <dbl>
```

3列目の食品群の列から群別の集計値を取り出せそうです。詳細な区分は今回は不要です。
6列目以降は栄養素別の列がずらっと並んでいるようです。今回はエネルギー、たんぱく質、脂質、炭水化物の４項目について調べようと思うので、その列だけ選択することにします。
これで必要な列の検討がつきました！

# 4.必要な行列だけを取り出す
まずは今回必要な項目を再度確認しましょ。
１：2017年のデータだけを扱う
２：グラフはエネルギー、たんぱく質、脂質、炭水化物の４つを作る
３：食品群別の集計値を使う（「穀類」は必要だけど、「米」や「米・加工品...」は不要）
４：ランキングを知りたいので、栄養素別に各食品群の割合を求める必要がある

です！
というわけで、まずは`filter()`と`select()`をつかって１と２に該当する箇所のデータだけ取り出すことにしましょう。
そして列名に【】が含まれているので、入力が面倒ですね。ついでに列名も変えちゃいましょう。

```r:行列の選択と列名変更
dat <- data %>% 
  filter(`時間軸(年次)` == "2017年") %>%  　　　　　　　　　　　　 # 2017年のデータだけ取り出す
  rename(                                　　　　　　　　　　　　# 列名変更　after = before
    food_gp = 食品群,
    energy = `エネルギー【kcal】`,
    protein_g = `たんぱく質【g】`,
    fat_g = `脂質【g】`,
    carbohydrate_g = `炭水化物【g】`
  ) %>% 
  select(food_gp, energy, protein_g, fat_g, carbohydrate_g)  # 必要な列だけ選択
```

```r:中身の確認
dat

> dat
# A tibble: 131 x 5
   food_gp                  energy protein_g fat_g carbohydrate_g
   <chr>                     <dbl>     <dbl> <dbl>          <dbl>
 1 総量                       1914      70.5  58.9          257. 
 2 動物性食品                  471      37.8  29             10.7
 3 植物性食品                 1443      32.7  29.9          246. 
 4 穀類                        760      15.2   4.7          158. 
 5 米・加工品                  513       7.7   1            113. 
 6 米                          506       7.5   0.9          112. 
 7 米加工品                      7       0.1   0              1.6
 8 小麦・加工品                229       7     3.6           41.2
 9 小麦粉類                     14       0.3   0.1            2.8
10 パン類（菓子パンを除く）     98       3.3   1.8           17.1
# … with 121 more rows
```
2017年のデータのうち、欲しい5列だけのデータセットができました！

続いて、食品群の集計値だけを取り出すことにします。
取り出したい要素はいくつもあるので、先にベクトルをつくってまとめちゃいます。

```r:食品群の集計値だけを取り出す
# 取り出したい要素のベクトルを作成
food_group <- c("総量", "穀類", "いも類", "砂糖・甘味料類", "豆類", "種実類",
                "野菜類", "果実類", "きのこ類", "藻類", "魚介類", "肉類",
                "卵類", "乳類", "油脂類", "菓子類", "嗜好飲料類", 
                "調味料・香辛料類", "その他")

# 取り出す
dat2 <- dat %>%
  filter(food_gp %in% food_group)
```

```r:中身チェック
dat2

> dat2
# A tibble: 18 x 5
   food_gp          energy protein_g fat_g carbohydrate_g
   <chr>             <dbl>     <dbl> <dbl>          <dbl>
 1 総量               1914      70.5  58.9          257. 
 2 穀類                760      15.2   4.7          158. 
 3 いも類               38       0.6   0.1            9.1
 4 砂糖・甘味料類       26       0     0              6.7
 5 豆類                 78       5.8   4.7            3.2
 6 種実類               15       0.5   1.3            0.7
 7 野菜類               74       3     0.5           16.7
 8 果実類               68       0.6   0.3           17.6
 9 きのこ類              3       0.4   0              1.1
10 藻類                  2       0.3   0              0.8
11 魚介類              107      13.1   4.9            1.7
12 肉類                204      15.3  14.8            0.4
13 卵類                 58       4.9   3.9            0.1
14 乳類                 92       4.5   4.4            8.4
15 油脂類              100       0    10.9            0  
16 菓子類               85       1.5   2.9           13.3
17 嗜好飲料類           90       1.1   0.1            8.1
18 調味料・香辛料類    112       3.8   5.5           11.1
```

# 5.栄養素別に各食品群の占める割合を求める

dat2は最初の行が総量になっています。この値を分母にして、各食品群の占める割合を求めます。

```r:各食品群の占める割合を求める
plot_data <- dat2 %>% 
  mutate(
    energy_per = (.$energy / .$energy[1] * 100),
    protein_per = (.$protein_g / .$protein_g[1] * 100),
    fat_per = (.$fat_g / .$fat_g[1] * 100),
    carbohydrate_per = (.$carbohydrate_g / .$carbohydrate_g[1] * 100)
  )
```

```r:中身チェック
plot_data

> plot_data
# A tibble: 18 x 9
   food_gp          energy protein_g fat_g carbohydrate_g energy_per protein_per fat_per carbohydrate_per
   <chr>             <dbl>     <dbl> <dbl>          <dbl>      <dbl>       <dbl>   <dbl>            <dbl>
 1 総量               1914      70.5  58.9          257.     100         100     100             100     
 2 穀類                760      15.2   4.7          158.      39.7        21.6     7.98           61.5   
 3 いも類               38       0.6   0.1            9.1      1.99        0.851   0.170           3.54  
 4 砂糖・甘味料類       26       0     0              6.7      1.36        0       0               2.61  
 5 豆類                 78       5.8   4.7            3.2      4.08        8.23    7.98            1.25  
 6 種実類               15       0.5   1.3            0.7      0.784       0.709   2.21            0.272 
 7 野菜類               74       3     0.5           16.7      3.87        4.26    0.849           6.50  
 8 果実類               68       0.6   0.3           17.6      3.55        0.851   0.509           6.85  
 9 きのこ類              3       0.4   0              1.1      0.157       0.567   0               0.428 
10 藻類                  2       0.3   0              0.8      0.104       0.426   0               0.311 
11 魚介類              107      13.1   4.9            1.7      5.59       18.6     8.32            0.662 
12 肉類                204      15.3  14.8            0.4     10.7        21.7    25.1             0.156 
13 卵類                 58       4.9   3.9            0.1      3.03        6.95    6.62            0.0389
14 乳類                 92       4.5   4.4            8.4      4.81        6.38    7.47            3.27  
15 油脂類              100       0    10.9            0        5.22        0      18.5             0     
16 菓子類               85       1.5   2.9           13.3      4.44        2.13    4.92            5.18  
17 嗜好飲料類           90       1.1   0.1            8.1      4.70        1.56    0.170           3.15  
18 調味料・香辛料類    112       3.8   5.5           11.1      5.85        5.39    9.34            4.32  
```

できたー！
このデータを元に棒グラフを作ってみます！

# 6.ggplotでグラフ化
長くなってきたので、ここはさくっといきます。とりあえずたんぱく質のグラフを作りましょう。
棒グラフは`geom_bar()`で作図します。
x軸を`食品群名(food_gp)`、y軸を`割合(protein_per)`にします。軸の値はx、yそれぞれの値を使いたいので`stat = "identity"`も忘れずに。

```r:たんぱく質のグラフ
plot_data %>% 
  ggplot(aes(food_gp, protein_per) +
  geom_bar(stat = "identity")
```
![graph1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1005648/c82dfa40-2988-c83a-7223-737905cbc98f.png)
なんか思ってたんと全然違うのができましたね。
x軸の値が文字化けして□になっちゃってるので何が何やら…てのと、並びがバラバラなので降順になって欲しいっすね。
`reorder()`で並びを指定しましょか。そして日本語のフォントも指定します。

```r:たんぱく質のグラフ（降順）
plot_data %>% 
  ggplot(aes(x = reorder(food_gp, -protein_per), y = protein_per)) +  # x軸はprotein_perに沿って降順
  geom_bar(stat = "identity") +
  theme_gray(base_family = "HiraKakuPro-W3")                          # 日本語フォントを指定
```
![graph2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1005648/b131a183-3108-f802-30b2-c0a0e4c0e937.png)
うーん、総量のバーが飛び出しててじゃまですね。そして軸タイトルもわかりやすくしたい。

```r:たんぱく質のグラフ（修正）
# 総量の行（1行目）を削除して上書き
plot_data <- plot_data[-1,]

plot_data %>% 
  ggplot(aes(x = reorder(food_gp, -protein_per), y = protein_per)) +
  geom_bar(stat = "identity") +
  theme_gray (base_family = "HiraKakuPro-W3") +
  labs(x = "食品群", y = "割合",　title = "たんぱく質摂取源となる食品群の摂取割合")  # タイトル追加・軸タイトル変更
```

![graph3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1005648/686bfe48-a8f4-ab78-862b-84efcad73cab.png)
ひとまず見やすくなりました！
やっぱりたんぱく質はお肉からとってますよね〜。
2番目は穀類ですってよ、さすが日本人。割合だけ見ると肉類とそんなに変わらないですね。
そして魚介類、豆類と続くようです。ふむふむ…。

他のグラフはコピペすればOKなので、本日はここまで。


# やってみた感想
おもしろかった！（雑）
他にもいろいろなことができそうですね〜。
APIからデータを引っ張ってくる方法も試してみたい。


# 参考元
- 政府統計の総合窓口(e-Stat)（https://www.e-stat.go.jp/）
- 国民健康・栄養調査（厚生労働省）(https://www.mhlw.go.jp/toukei/itiran/gaiyo/k-eisei.html)
- R Programming Tutorial - Learn the Basics of Statistical Computing (https://youtu.be/_V8eKsto3Ug)

`pacman`のパッケージは上記の動画で覚えました