7  オープンデータの整形

オープンデータとは、許可される範囲内で誰もが自由に利用・再利用・再配布できるデータのことを指す。日本でも、政府や地方行政、研究機関が統計調査の結果を公開する例が増えているが、公開されたデータはRでの処理をすぐに実行できない状態であることがしばしばある。これらのファイルもtidyverseのパッケージを使うことで読み込みからデータ整形、可視化までの作業を一貫して行うことが可能である。本章では、実際に公開されているオープンデータを例にtidyverseによる作業の例を紹介しよう。

# 必要なパッケージを準備する
library(tidyverse)

── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

✔ ggplot2 3.3.5     ✔ purrr   0.3.4
✔ tibble  3.1.6     ✔ dplyr   1.0.9
✔ tidyr   1.2.0     ✔ stringr 1.4.0
✔ readr   2.1.2     ✔ forcats 0.5.1

── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

library(lubridate)

Attaching package: 'lubridate'

The following objects are masked from 'package:base':

    date, intersect, setdiff, union

7.1 政府統計の総合窓口 e-Stat

政府統計の総合窓口 (e-Stat) は日本国内の政府統計のポータルサイトである。各府省等が公表する各種の統計情報を閲覧、ダウンロード可能であり、Rなどのプログラムを介したデータ参照を行うためのAPIも整備されている。

Rからはe-StatのAPIを操作するestatapiパッケージを使うことで、提供されるデータの一覧や、データの読み込み、ファイルダウンロードが可能となる。APIの利用にはe-Statでのメールアドレス等の登録が必要となるが、ここではあらかじめ用意したアプリケーションIDを用いる例を紹介する。なお本書の付録として、e-Statからダウンロードしたデータを掲載しているので、アプリケーションIDの取得を行わないでも実行可能である。

library(estatapi)

このサービスは、政府統計総合窓口(e-Stat)のAPI機能を使用していますが、サービスの内容は国によって保証されたものではありません。

7.1.1 学校保健統計調査

「学校保健統計調査」は、学校における幼児、児童及び生徒の発育及び健康の状態を明らかにすることを目的に行われている。調査の対象は、文部科学大臣があらかじめ指定する学校に在籍する満5歳から17歳（4月1日現在）までの幼児、児童及び生徒となっている。調査は毎年実施され、学校保健安全法により義務づけられている健康診断の結果に基づき、発育及び健康状態に関する事項（身長、体重及び被患率等）に関する調査結果が記録されている。

このデータがe-StatのAPIで取得可能かどうか、まずはestatapiパッケージのestat_getStatsList()で調べてみよう。estat_getStatsList()には引数appIdがあり、ここに取得したアプリケーションIDを指定する。

df_list <- 
  estat_getStatsList(appId      = "<APIで取得したアプリケーションID>", 
                     searchWord = "学校保健統計調査")

df_listを見ると、学校保健統計調査に関する様々なデータがあることが確認できる。その中から、今回は「学校保健統計調査による身体発育値」の名目で記録されたデータを抽出してみよう。

df_list %>% 
  filter(str_detect(TITLE, regex("学校保健統計調査による身体発育値")))

# A tibble: 5 × 13
  `@id`      STAT_NAME GOV_ORG STATISTICS_NAME TITLE CYCLE SURVEY_DATE OPEN_DATE
  <chr>      <chr>     <chr>   <chr>           <chr> <chr> <chr>       <chr>    
1 0003079318 学校保健… 文部科… 学校保健統計調… （参… 年度… 0           2014-03-…
2 0003107473 学校保健… 文部科… 学校保健統計調… （参… 年度… 201404-201… 2015-03-…
3 0003107409 学校保健… 文部科… 学校保健統計調… （参… 年度… 201404-201… 2015-03-…
4 0003146920 学校保健… 文部科… 学校保健統計調… （参… 年度… 201504-201… 2016-03-…
5 0003146941 学校保健… 文部科… 学校保健統計調… （参… 年度… 201504-201… 2016-03-…
# … with 5 more variables: SMALL_AREA <chr>, MAIN_CATEGORY <chr>,
#   SUB_CATEGORY <chr>, OVERALL_TOTAL_NUMBER <chr>, UPDATED_DATE <chr>

「学校保健統計調査による身体発育値」の項目には5件のデータが該当している。STATISTICS_NAME列を見ると、データは平成22年度以降、平成26年度、平成27年度の3期間に渡り、平成26年と27年ではLMS法を採用した項目が別にあることが確認できる。

estatapiパッケージを使ったe-Statのデータ取得には、対象データの第一列に記録された”@id”の値が必要となる。そのため、取得対象のデータを決めたらこの値を控えておくと良い。次の処理では、データフレーム中の”@id”列の値を参照するが、LMS法の項目を区別するために正規表現による抽出を試みている。

# 「学校保健統計調査による身体発育値」のデータを取得するためのIDを控えておく
ids <- 
  df_list %>% 
  # regex()内で正規表現のパターンマッチを行い、「身体発育値」で終わる項目を抽出する
  filter(str_detect(TITLE, regex("学校保健統計調査による身体発育値$"))) %>% 
  .$`@id`

ids

[1] "0003079318" "0003107473" "0003146920"

抽出されたデータの件数は3件である。次はこの”@id”から、R上にデータを読み込んでみよう。データ取得はestat_getStatsData()で対象の統計データIDおよびアプリケーションIDを指定して行う。

estat_getStatsData(
      appId = config::get("estat_token"),
      statsDataId = ids[1])

学校保健統計調査は年に一度であるが、今回取得可能なデータの件数は3期間分ある。これを一つのデータフレームとして処理するには、次のようにpurrr::map_dfr()を使い、estat_getStatsData()のstatsDataId引数にidsの要素を一つずつ渡して実行し、行方向への結合を行うと良い。

df_age_height <- 
  ids %>% 
  map_dfr(
    ~ estat_getStatsData(
      appId = "<API token>",
      statsDataId = .x
    )
  )

それではtibble::glimpse()を使い取得したデータを確認しよう。

glimpse(df_age_height)

Rows: 2,184
Columns: 14
$ tab_code                     <chr> "0000300", "0000300", "0000300", "0000300…
$ 表章項目                     <chr> "身長のパーセンタイル値(cm)", "身長のパー…
$ cat01_code                   <chr> "20", "20", "20", "20", "20", "20", "20",…
$ 性別                         <chr> "男", "男", "男", "男", "男", "男", "男",…
$ cat03_code                   <chr> "00000100", "00000100", "00000100", "0000…
$ パーセンタイル               <chr> "3%", "3%", "3%", "3%", "3%", "3%", "3%",…
$ cat06_code                   <chr> "10", "10", "20", "20", "30", "30", "40",…
$ `学校種別-年齢別（5～17歳）` <chr> "幼稚園（5歳）", "幼稚園（5歳）", "小学校…
$ area_code                    <chr> "00000", "00000", "00000", "00000", "0000…
$ 都道府県別                   <chr> "全国", "全国", "全国", "全国", "全国", "…
$ time_code                    <chr> "2013000000", "2012000000", "2013000000",…
$ `時間軸(年次)`               <chr> "2013年", "2012年", "2013年", "2012年", "…
$ value                        <dbl> 101.7, 101.6, 107.5, 107.7, 113.0, 113.3,…
$ `時間軸（年度次）`           <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…

次に、ここからdplyrのデータ操作関数を使いながらデータを加工していこう。まずは今後使う列だけを選択する。また、この際に日本語の列名から適当な列名に変更を行う。各列の値はそれぞれ以下の項目を示している。

type: 計測項目。身長、体重、座高の3種がある。

gender: 性別。「男」と「女」の二値が記録されている。

percent: データのパーセンタイル値。

age: 対象者の学年および年齢。

fiscal_year: 調査の実施年度。

value: 計測項目 typeごとの計測値。

df_age_height <- 
  df_age_height %>% 
  select(type = `表章項目`,
         gender = `性別`,
         percent = `パーセンタイル`,
         age = `学校種別-年齢別（5～17歳）`,
         fiscal_year = `時間軸(年次)`,
         value)

glimpse(df_age_height)

Rows: 2,184
Columns: 6
$ type        <chr> "身長のパーセンタイル値(cm)", "身長のパーセンタイル値(cm)"…
$ gender      <chr> "男", "男", "男", "男", "男", "男", "男", "男", "男", "男"…
$ percent     <chr> "3%", "3%", "3%", "3%", "3%", "3%", "3%", "3%", "3%", "3%"…
$ age         <chr> "幼稚園（5歳）", "幼稚園（5歳）", "小学校（6歳）", "小学校…
$ fiscal_year <chr> "2013年", "2012年", "2013年", "2012年", "2013年", "2012年"…
$ value       <dbl> 101.7, 101.6, 107.5, 107.7, 113.0, 113.3, 118.2, 118.6, 12…

type、ageの列には、括弧の中に計測値の単位、対象者の年齢が記録されている。これをtidyrのextract()で新たな列として独立させよう。

df_age_height <- 
  df_age_height %>% 
  extract(type, into = c("type", "unit"),
           regex = "(.+)(\\([[:alnum:]]{1,}+\\))") %>% 
  extract(age, into = c("school", "age"),
          regex = "(^.+)(（.+$)")

glimpse(df_age_height)

Rows: 2,184
Columns: 8
$ type        <chr> "身長のパーセンタイル値", "身長のパーセンタイル値", "身長…
$ unit        <chr> "(cm)", "(cm)", "(cm)", "(cm)", "(cm)", "(cm)", "(cm)", "(…
$ gender      <chr> "男", "男", "男", "男", "男", "男", "男", "男", "男", "男"…
$ percent     <chr> "3%", "3%", "3%", "3%", "3%", "3%", "3%", "3%", "3%", "3%"…
$ school      <chr> "幼稚園", "幼稚園", "小学校", "小学校", "小学校", "小学校"…
$ age         <chr> "（5歳）", "（5歳）", "（6歳）", "（6歳）", "（7歳）", "（…
$ fiscal_year <chr> "2013年", "2012年", "2013年", "2012年", "2013年", "2012年"…
$ value       <dbl> 101.7, 101.6, 107.5, 107.7, 113.0, 113.3, 118.2, 118.6, 12…

次に、単位のついたageとfiscal_yearの列を対象として、数値を抽出する処理を実行する。ここでは対象列の指定にmutate_at()、数値の抽出にparse_number()を適用する。また、mutate_at()の適用後も、変数は元の文字列のデータ型を保持しているため、これをtype_convert()を使い数値変数として扱えるようにしておこう。

df_age_height <- 
  df_age_height %>% 
  mutate_at(vars(c("age", "fiscal_year")), parse_number) %>% 
  type_convert()

── Column specification ────────────────────────────────────────────────────────
cols(
  type = col_character(),
  unit = col_character(),
  gender = col_character(),
  percent = col_character(),
  school = col_character()
)

df_age_height

# A tibble: 2,184 × 8
   type                   unit  gender percent school   age fiscal_year value
   <chr>                  <chr> <chr>  <chr>   <chr>  <dbl>       <dbl> <dbl>
 1 身長のパーセンタイル値 (cm)  男     3%      幼稚園     5        2013  102.
 2 身長のパーセンタイル値 (cm)  男     3%      幼稚園     5        2012  102.
 3 身長のパーセンタイル値 (cm)  男     3%      小学校     6        2013  108.
 4 身長のパーセンタイル値 (cm)  男     3%      小学校     6        2012  108.
 5 身長のパーセンタイル値 (cm)  男     3%      小学校     7        2013  113 
 6 身長のパーセンタイル値 (cm)  男     3%      小学校     7        2012  113.
 7 身長のパーセンタイル値 (cm)  男     3%      小学校     8        2013  118.
 8 身長のパーセンタイル値 (cm)  男     3%      小学校     8        2012  119.
 9 身長のパーセンタイル値 (cm)  男     3%      小学校     9        2013  123.
10 身長のパーセンタイル値 (cm)  男     3%      小学校     9        2012  123.
# … with 2,174 more rows

こうして加工した学校保健統計調査のデータを、最後にggplot2を使って可視化してみよう。まずは3種の計測項目(“type”)のデータのばらつきをgeom_density()で表示させる。複数の項目からなるデータを、項目別に描画するには次のようにfacet_wrap()を使うと良い。

df_age_height %>% 
  ggplot(aes(value, color = gender)) + 
  geom_density() +
  # typeごとに図を分ける。
  facet_wrap(~ type, 
             # x軸はtypeごとに独立した軸をもつことを指定する
             scales = "free_x",
             # 分割して描画する図の列数
             ncol = 3)

次の図は、2013年に記録された男女の身長データについてプロットしたものである。まずfilter()を使い対象のデータを抽出する。x軸には年齢を指定するが、ここで元の”age”が数値であるためにforcats::as_factor()を使い因子へと変換させた。

df_age_height %>%
  filter(str_detect(type, "身長"), 
         fiscal_year == 2013) %>% 
  ggplot(aes(forcats::as_factor(as.character(age)), value, 
             group = percent, color = percent)) +
  geom_point() +
  geom_line() +
  xlab("年齢") +
  ylab("身長 (cm)")

  facet_wrap(~ gender, ncol = 1)

<ggproto object: Class FacetWrap, Facet, gg>
    compute_layout: function
    draw_back: function
    draw_front: function
    draw_labels: function
    draw_panels: function
    finish_data: function
    init_scales: function
    map_data: function
    params: list
    setup_data: function
    setup_params: function
    shrink: TRUE
    train_scales: function
    vars: function
    super:  <ggproto object: Class FacetWrap, Facet, gg>

7.2 気象庁の過去データ

気象庁のウェブサイト ()では、過去の起床データをダウンロードできるサービスが提供されている。ここでは、観測地点「つくば（館野）」での2017年における気象データを用意した。このデータを使い、気象庁のデータ整形と簡単な可視化を行おう。このデータを使わなくても気象庁の「過去の気象データ・ダウンロード」ページでダウンロードしたファイルには同様の処理が有効だろう。

まずダウンロードしたファイルの読み込みを実施する。元のデータはcsv形式で提供されるが、エンコーディングがWindows向けとなっているため、文字化けを防ぐためにreadr::locale()でエンコードの指定する。また先頭行はファイルをダウンロードした時刻と対象の観測所名が記録されているので、これらの読み込みを引数skipで除外しよう。

df_weather <- 
  read_csv("data-raw/jma_tsukubka_2017data.csv", 
           locale = locale(encoding = "cp932"), 
           skip = 3)

New names:
Rows: 367 Columns: 24
── Column specification
──────────────────────────────────────────────────────── Delimiter: "," chr
(17): 年月日, 平均気温(℃)...3, 平均気温(℃)...4, 最低気温(℃)...6,
最低気温(℃)...7, 最高気温(℃)..... dbl (6): 平均気温(℃)...2, 最低気温(℃)...5,
最高気温(℃)...8, 降水量の合計(mm)...11, 日照時間(時間)..... lgl (1):
平均雲量(10分比)...22
ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
Specify the column types or set `show_col_types = FALSE` to quiet this message.
• `平均気温(℃)` -> `平均気温(℃)...2`
• `平均気温(℃)` -> `平均気温(℃)...3`
• `平均気温(℃)` -> `平均気温(℃)...4`
• `最低気温(℃)` -> `最低気温(℃)...5`
• `最低気温(℃)` -> `最低気温(℃)...6`
• `最低気温(℃)` -> `最低気温(℃)...7`
• `最高気温(℃)` -> `最高気温(℃)...8`
• `最高気温(℃)` -> `最高気温(℃)...9`
• `最高気温(℃)` -> `最高気温(℃)...10`
• `降水量の合計(mm)` -> `降水量の合計(mm)...11`
• `降水量の合計(mm)` -> `降水量の合計(mm)...12`
• `降水量の合計(mm)` -> `降水量の合計(mm)...13`
• `降水量の合計(mm)` -> `降水量の合計(mm)...14`
• `日照時間(時間)` -> `日照時間(時間)...15`
• `日照時間(時間)` -> `日照時間(時間)...16`
• `日照時間(時間)` -> `日照時間(時間)...17`
• `日照時間(時間)` -> `日照時間(時間)...18`
• `平均風速(m/s)` -> `平均風速(m/s)...19`
• `平均風速(m/s)` -> `平均風速(m/s)...20`
• `平均風速(m/s)` -> `平均風速(m/s)...21`
• `平均雲量(10分比)` -> `平均雲量(10分比)...22`
• `平均雲量(10分比)` -> `平均雲量(10分比)...23`
• `平均雲量(10分比)` -> `平均雲量(10分比)...24`

データを確認すると、平均気温、最高気温などの各項目にその値と品質情報、均質番号が記録された列があることがわかる。欲しいのは項目の値であるので品質情報および均質番号の列はデータから削除しよう。これには、品質情報および均質番号の列が最初の行だけ値があり、あとは欠損値であることを利用して、対象の変数をベクトルとして得る次の処理を実行すると良いだろう。

glimpse(df_weather)

Rows: 367
Columns: 24
$ 年月日                  <chr> NA, NA, "2017/1/1", "2017/1/2", "2017/1/3", "2…
$ `平均気温(℃)...2`       <dbl> NA, NA, 4.0, 3.9, 4.8, 5.7, 3.6, 0.7, 0.8, 3.0…
$ `平均気温(℃)...3`       <chr> NA, "品質情報", "8", "8", "8", "8", "8", "8", …
$ `平均気温(℃)...4`       <chr> NA, "均質番号", "1", "1", "1", "1", "1", "1", …
$ `最低気温(℃)...5`       <dbl> NA, NA, -3.0, -1.4, -2.2, -2.4, -3.2, -5.6, -5…
$ `最低気温(℃)...6`       <chr> NA, "品質情報", "8", "8", "8", "8", "8", "8", …
$ `最低気温(℃)...7`       <chr> NA, "均質番号", "1", "1", "1", "1", "1", "1", …
$ `最高気温(℃)...8`       <dbl> NA, NA, 13.3, 12.5, 13.8, 14.1, 9.7, 8.3, 8.8,…
$ `最高気温(℃)...9`       <chr> NA, "品質情報", "8", "8", "8", "8", "8", "8", …
$ `最高気温(℃)...10`      <chr> NA, "均質番号", "1", "1", "1", "1", "1", "1", …
$ `降水量の合計(mm)...11` <dbl> NA, NA, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 9.5…
$ `降水量の合計(mm)...12` <chr> NA, "現象なし情報", "1", "1", "1", "1", "1", "…
$ `降水量の合計(mm)...13` <chr> NA, "品質情報", "8", "8", "8", "8", "8", "8", …
$ `降水量の合計(mm)...14` <chr> NA, "均質番号", "1", "1", "1", "1", "1", "1", …
$ `日照時間(時間)...15`   <dbl> NA, NA, 9.5, 4.9, 9.3, 9.2, 8.9, 8.1, 9.3, 0.3…
$ `日照時間(時間)...16`   <chr> NA, "現象なし情報", "0", "0", "0", "0", "0", "…
$ `日照時間(時間)...17`   <chr> NA, "品質情報", "8", "8", "8", "8", "8", "8", …
$ `日照時間(時間)...18`   <chr> NA, "均質番号", "1", "1", "1", "1", "1", "1", …
$ `平均風速(m/s)...19`    <dbl> NA, NA, 1.4, 1.3, 2.0, 2.3, 1.8, 1.3, 1.1, 2.2…
$ `平均風速(m/s)...20`    <chr> NA, "品質情報", "8", "8", "8", "8", "8", "8", …
$ `平均風速(m/s)...21`    <chr> NA, "均質番号", "1", "1", "1", "1", "1", "1", …
$ `平均雲量(10分比)...22` <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `平均雲量(10分比)...23` <chr> NA, "品質情報", "0", "0", "0", "0", "0", "0", …
$ `平均雲量(10分比)...24` <chr> NA, "均質番号", "1", "1", "1", "1", "1", "1", …

select_vars <-
  df_weather %>%
  slice(2L) %>%
  select_if(.predicate = is.na(.)) %>%
  names()

df_weather <- 
  df_weather %>% 
  slice(-c(1:2L)) %>% 
  select(select_vars)

Note: Using an external vector in selections is ambiguous.
ℹ Use `all_of(select_vars)` instead of `select_vars` to silence this message.
ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
This message is displayed once per session.

df_weather <- 
  df_weather %>% 
  type_convert() %>% 
  rename(date = `年月日`) %>% 
  mutate(date = ymd(date))

── Column specification ────────────────────────────────────────────────────────
cols(
  年月日 = col_character()
)

df_weather %>% 
  # extract(`平均気温(℃)`, into = "temp_mean", convert = TRUE) %>% 
  # extract(`最低気温(℃)`, into = "temp_min", convert = TRUE) %>% 
  # extract(`最高気温(℃)`, into = "temp_max", convert = TRUE) %>% 
  gather(type, value, -date, convert = TRUE) %>% 
  type_convert()

── Column specification ────────────────────────────────────────────────────────
cols(
  type = col_character()
)

# A tibble: 2,555 × 3
   date       type            value
   <date>     <chr>           <dbl>
 1 2017-01-01 平均気温(℃)...2   4  
 2 2017-01-02 平均気温(℃)...2   3.9
 3 2017-01-03 平均気温(℃)...2   4.8
 4 2017-01-04 平均気温(℃)...2   5.7
 5 2017-01-05 平均気温(℃)...2   3.6
 6 2017-01-06 平均気温(℃)...2   0.7
 7 2017-01-07 平均気温(℃)...2   0.8
 8 2017-01-08 平均気温(℃)...2   3  
 9 2017-01-09 平均気温(℃)...2   5.6
10 2017-01-10 平均気温(℃)...2   6  
# … with 2,545 more rows

df_weather_temp <- 
  df_weather %>% 
  select(1, contains("気温"))

df_weather_temp %>% 
  gather(class, value, -date) %>% 
  ggplot(aes(date, value, color = class)) +
  geom_line()

7.3 「青果物卸売市場調査（平成29年年間計及び月別結果）」（農林水産省）

library(readxl)

xx <- c("卸売数量", "卸売価額", "卸売価格")

df_1 <- 
  read_excel("data-raw/seika_oroshi17.xls", sheet = 7, skip = 5, 
             col_names = c(NA, 
                           "国産_輸入", "品目", "種別",
                           NA,
                           str_c("主要都市の市場計", "_", xx),
                           str_c("対　　前　　年　　比", "_", xx),
                           str_c("１月", "_", xx),
                           str_c("２月", "_", xx),
                           str_c("３月", "_", xx),
                           NA)) %>% 
  mutate_at(vars(contains("月")), .funs = funs(recode(., `-` = NA_character_))) %>% 
  slice(-c(1:4)) %>% 
  select(-num_range("X__", 1:3)) %>% 
  fill(`国産_輸入`, `品目`, .direction = "down") %>% 
  drop_na(`国産_輸入`)

df_2 <- 
  read_excel("data-raw/seika_oroshi17.xls", sheet = 8, skip = 5,
             col_names = c(NA, 
                           "国産_輸入", "品目", "種別",
                           NA,
                           str_c("４月", "_", xx),
                           str_c("５月", "_", xx),
                           str_c("６月", "_", xx),
                           str_c("７月", "_", xx),
                           str_c("８月", "_", xx),
                           str_c("９月", "_", xx),
                           NA)) %>% 
  mutate_at(vars(contains("月")), .funs = funs(recode(., `-` = NA_character_))) %>% 
  slice(-c(1:4)) %>% 
  select(-num_range("X__", 1:3)) %>% 
  fill(`国産_輸入`, `品目`, .direction = "down") %>% 
  drop_na(`国産_輸入`)

df_3 <- 
  read_excel("data-raw/seika_oroshi17.xls", sheet = 9, skip = 5,
           col_names = c(NA, "国産_輸入", "品目", "種別",
                         NA,
                         str_c("10月", "_", xx),
                         str_c("11月", "_", xx),
                         str_c("12月", "_", xx))) %>% 
  mutate_at(vars(contains("月")), .funs = funs(recode(., `-` = NA_character_))) %>% 
  slice(-c(1:4)) %>% 
  select(-num_range("X__", 1:2)) %>% 
  fill(`国産_輸入`, `品目`, .direction = "down") %>% 
  drop_na(`国産_輸入`)

df_seika <- 
  left_join(df_1, df_2, by = c("国産_輸入", "品目", "種別")) %>% 
  left_join(df_3, by = c("国産_輸入", "品目", "種別")) %>% 
  filter_at(vars(`品目`, `種別`), any_vars(!is.na(.))) %>% 
  group_by_at(vars(-contains("主要都市の市場計"), -contains("対　　前　　年　　比"))) %>% 
  nest() %>%
  gather(key, value, -`国産_輸入`, -`品目`, -`種別`, -data) %>% 
  separate(key, c("month", "key"), sep = "_") %>% 
  type_convert()

df_seika %>% 
  filter(key == "卸売数量", month == "１月") %>% 
  group_by(`品目`) %>% 
  summarise(value = sum(value)) %>% 
  mutate(`品目` = forcats::fct_reorder(`品目`, value)) %>% 
  drop_na() %>% 
  ggplot(aes(`品目`, value)) +
  geom_bar(stat = "identity")

df_seika %>% 
  filter(key == "卸売数量", `品目` == "りんご計", !is.na(種別)) %>% 
  ggplot(aes(month, value, color = `種別`, group = `種別`)) +
  geom_point() +
  geom_line() +
  facet_wrap(~ `国産_輸入`)