だんだん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/
キーワード検索から「国民健康・栄養調査」を検索します。
栄養素ごとに、どの食品からどれくらいの量を食べているのかがわかるデータ（食品群別栄養素等摂取量）を手に入れようと思います。
検索結果はひとつしか出てこないので、表示されるままクリックしてここまできます。
データベースをクリックして 年次 → 2017年 → 表番号９のDB（画像赤枠） の順に進みます。
この画面が表示されたら、特になにもいじらずそのまま右のダウンロードをクリック。
表ダウンロードのページに飛ぶので、ヘッダの出力：出力しない、コードの出力：出力しない を選んで（他はそのまま）、下部に表示されるダウンロードボタンを押します。
お好みの場所、ファイル名で保存します。
データがDLできたら、ファイルを開いて中身を確認してみます。
Macだと勝手にNumbersで開いちゃうけど、ここはExcelで開きます
ちゃんとA1からデータで埋まっています。やったぜ
これで必要なデータが入手できました！
2.DLしたデータの読み込み
DLしたデータをRに読み込んでいきましょう。まずはtidyverseとmagrittrを召喚。
library(tidyverse)
library(magrittr)
わたしはlibraryを2回呼び出すのが面倒なので、pacmanというパッケージで一気に召喚してます。
# pacmanが入ってなかったらインストール
if(!(require(pacman))) install.packages("pacman")
# tidyverseとmagrittrを呼び出す
pacman::p_load(pacman, tidyverse, magrittr)
readr()でDLしたデータを読み込むと…
data <- read_r('ファイル名.csv')
> data <- read_csv('ファイル名.csv')
make.names(x) でエラー:
'<94>N<97><ee><8a><4b>��' に不正なマルチバイト文字があります
ってな感じでエラーが出ます。
不正なマルチバイト文字ってことなので、文字のエンコードを指定する必要がありそうです。
日本語のデータなので、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出力されているはずです。
ヘッダーを出力していると、読み込みの際にエラーが出るので1行目から11行目をスキップする必要があります。
skip = を追加してどこまでスキップする行数を指定すればOKです。
# 11行目まではスキップして12行目から読み込んでね
data <- read_csv('ファイル名.csv', locale = locale(encoding = "SHIFT-JIS"), skip = 11)
3.読み込んだデータをざくっと眺める
読み込んだデータの中身をstr()で確認してみます。
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()を使ってさらに確認してみます。
data$年齢階級 %>% summary.factor()
> data$年齢階級 %>% summary.factor()
20歳以上
262
年齢階級の列は同じ情報だけのようなので特に気にする必要はなさそうです。
data$`時間軸(年次)` %>% summary.factor()
> data$`時間軸(年次)` %>% summary.factor()
2016年 2017年
131 131
あれ？2017年のデータをDLしたつもりが、2016年の情報も入ってるみたいですね。
今回はとりあえず、2017年のデータだけを使いましょう。
slice()を使って最初の数行の内容も確認してみます。
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()をつかって１と２に該当する箇所のデータだけ取り出すことにしましょう。
そして列名に【】が含まれているので、入力が面倒ですね。ついでに列名も変えちゃいましょう。
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) # 必要な列だけ選択
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列だけのデータセットができました！
続いて、食品群の集計値だけを取り出すことにします。
取り出したい要素はいくつもあるので、先にベクトルをつくってまとめちゃいます。
# 取り出したい要素のベクトルを作成
food_group <- c("総量", "穀類", "いも類", "砂糖・甘味料類", "豆類", "種実類",
"野菜類", "果実類", "きのこ類", "藻類", "魚介類", "肉類",
"卵類", "乳類", "油脂類", "菓子類", "嗜好飲料類",
"調味料・香辛料類", "その他")
# 取り出す
dat2 <- dat %>%
filter(food_gp %in% food_group)
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は最初の行が総量になっています。この値を分母にして、各食品群の占める割合を求めます。
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)
)
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"も忘れずに。
plot_data %>%
ggplot(aes(food_gp, protein_per) +
geom_bar(stat = "identity")
なんか思ってたんと全然違うのができましたね。
x軸の値が文字化けして□になっちゃってるので何が何やら…てのと、並びがバラバラなので降順になって欲しいっすね。
reorder()で並びを指定しましょか。そして日本語のフォントも指定します。
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") # 日本語フォントを指定
うーん、総量のバーが飛び出しててじゃまですね。そして軸タイトルもわかりやすくしたい。
# 総量の行（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 = "たんぱく質摂取源となる食品群の摂取割合") # タイトル追加・軸タイトル変更
ひとまず見やすくなりました！
やっぱりたんぱく質はお肉からとってますよね〜。
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のパッケージは上記の動画で覚えました