Lesson 3-1: dplyr入門

2018/11/04

前のステップ

Lesson 2-2: データフレーム・tibble では次のことを勉強しました。

このセッションの目標

Tidyverse

このセッションでは，tidyverse というパッケージ群を使います。

library(tidyverse)
## ── Attaching packages ───────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.0.0     ✔ purrr   0.2.5
## ✔ tibble  1.4.2     ✔ dplyr   0.7.6
## ✔ tidyr   0.8.1     ✔ stringr 1.3.1
## ✔ readr   1.1.1     ✔ forcats 0.3.0
## ── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

メッセージを一度読んでおきましょう。library(tidyverse) というコマンドは（今後増減するかもしれませんが）現在 8個のパッケージを読み込みます。

本日のメインは dplyr （ディープライアーと読む）です。

“Conflicts” という少し不穏なメッセージも出ていますが，これはすでに読み込まれている関数が上書きされたことに対する警告です。

dplyr::filter() が stats::filter() を上書きしたので，パッケージ名を付けずに filter() を呼び出すと dplyr::filter() が 呼ばれます。もし，stats::filter() を使いたい場合には stats:: を付けます。lag() についても同様です。

Penn World Table

本日は，Penn World Table という主要なマクロ指標を国際比較できるデータベースを使います。1950〜2014年，182カ国のデータが掲載されています。 データの詳細については以下の論文を参照してください。

Feenstra, Robert C., Robert Inklaar and Marcel P. Timmer (2015), “The Next Generation of the Penn World Table” American Economic Review, 105(10), 3150-3182, available for download at http://www.ggdc.net/pwt

ダウンロード

まずは本日使うデータをダウンロードしましょう。いつも通り，

  • RClub 用のプロジェクトを開いている状態で，
  • その中に Data フォルダがあること

を想定しています。

次のコマンドをコンソールで実行してください。（コピー＆ペーストしてください）

download.file(url = "http://www.rug.nl/ggdc/docs/pwt90.dta", 
              destfile = "Data/pwt90.dta", mode = "wb")

これは Penn World Table version 9.0 の Stata 形式のデータをダウンロードします。

読み込み

読み込むときは

pwt90 <- haven::read_dta("Data/pwt90.dta")

とします。haven というSTATA や SPSS のデータ形式を読み込むためのパッケージを使っています。読み込んだデータに pwt90 という名前を付けたので，コンソールに pwt90 と打ち込むと最初の方にある行と列が表示されます。

pwt90
## # A tibble: 11,830 x 47
##    countrycode country currency_unit  year rgdpe rgdpo   pop   emp   avh
##    <chr>       <chr>   <chr>         <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1 ABW         Aruba   Aruban Guild…  1950    NA    NA    NA    NA    NA
##  2 ABW         Aruba   Aruban Guild…  1951    NA    NA    NA    NA    NA
##  3 ABW         Aruba   Aruban Guild…  1952    NA    NA    NA    NA    NA
##  4 ABW         Aruba   Aruban Guild…  1953    NA    NA    NA    NA    NA
##  5 ABW         Aruba   Aruban Guild…  1954    NA    NA    NA    NA    NA
##  6 ABW         Aruba   Aruban Guild…  1955    NA    NA    NA    NA    NA
##  7 ABW         Aruba   Aruban Guild…  1956    NA    NA    NA    NA    NA
##  8 ABW         Aruba   Aruban Guild…  1957    NA    NA    NA    NA    NA
##  9 ABW         Aruba   Aruban Guild…  1958    NA    NA    NA    NA    NA
## 10 ABW         Aruba   Aruban Guild…  1959    NA    NA    NA    NA    NA
## # ... with 11,820 more rows, and 38 more variables: hc <dbl>, ccon <dbl>,
## #   cda <dbl>, cgdpe <dbl>, cgdpo <dbl>, ck <dbl>, ctfp <dbl>,
## #   cwtfp <dbl>, rgdpna <dbl>, rconna <dbl>, rdana <dbl>, rkna <dbl>,
## #   rtfpna <dbl>, rwtfpna <dbl>, labsh <dbl>, delta <dbl>, xr <dbl>,
## #   pl_con <dbl>, pl_da <dbl>, pl_gdpo <dbl>, i_cig <dbl+lbl>,
## #   i_xm <dbl+lbl>, i_xr <dbl+lbl>, i_outlier <dbl+lbl>, cor_exp <dbl>,
## #   statcap <dbl>, csh_c <dbl>, csh_i <dbl>, csh_g <dbl>, csh_x <dbl>,
## #   csh_m <dbl>, csh_r <dbl>, pl_c <dbl>, pl_i <dbl>, pl_g <dbl>,
## #   pl_x <dbl>, pl_m <dbl>, pl_k <dbl>

11,830行，47列あることも分かります。列は変数に対応しています。

STATA のラベル

pwt90.dta はもともとSTATAのデータなので，変数には説明用のラベルがついています。haven::read_dta() はラベル情報も同時にインポートします。単独の変数（列）についたラベルを知りたい場合は，次のコマンドを実行してみてください。

attr(pwt90$rgdpo, "label")
## [1] "Output-side real GDP at chained PPPs (in mil. 2011US$)"

すべての変数とラベルの関係を一覧するコードブックを作っておくと便利かもしれません。例えば，次のようにします。（map_chr() は purrr の関数で，同じ操作をベクトルの各要素やデータフレームの各列に繰り返すときに使います。何をやっているかを今すぐ理解できなくても構いません。）

labs <- map_chr(colnames(pwt90), ~ attr(pwt90[[.x]], "label"))  # ラベルの取得
(codebook <- tibble(variable = colnames(pwt90), label = labs))  # tibble に変換
## # A tibble: 47 x 2
##    variable     label                                                     
##    <chr>        <chr>                                                     
##  1 countrycode  3-letter ISO country code                                 
##  2 country      Country name                                              
##  3 currency_un… Currency unit                                             
##  4 year         Year                                                      
##  5 rgdpe        Expenditure-side real GDP at chained PPPs (in mil. 2011US…
##  6 rgdpo        Output-side real GDP at chained PPPs (in mil. 2011US$)    
##  7 pop          Population (in millions)                                  
##  8 emp          Number of persons engaged (in millions)                   
##  9 avh          Average annual hours worked by persons engaged (source: T…
## 10 hc           Human capital index, see note hc                          
## # ... with 37 more rows

全体を見るためには View() が使えます。

View(codebook)
variable label
countrycode 3-letter ISO country code
country Country name
currency_unit Currency unit
year Year
rgdpe Expenditure-side real GDP at chained PPPs (in mil. 2011US$)
rgdpo Output-side real GDP at chained PPPs (in mil. 2011US$)
pop Population (in millions)
emp Number of persons engaged (in millions)
avh Average annual hours worked by persons engaged (source: The Conference Board)
hc Human capital index, see note hc
ccon Real consumption of households and government, at current PPPs (in mil. 2011US$)
cda Real domestic absorption, see note cda
cgdpe Expenditure-side real GDP at current PPPs (in mil. 2011US$)
cgdpo Output-side real GDP at current PPPs (in mil. 2011US$)
ck Capital stock at current PPPs (in mil. 2011US$)
ctfp TFP level at current PPPs (USA=1)
cwtfp Welfare-relevant TFP levels at current PPPs (USA=1)
rgdpna Real GDP at constant 2011 national prices (in mil. 2011US$)
rconna Real consumption at constant 2011 national prices (in mil. 2011US$)
rdana Real domestic absorption at constant 2011 national prices (in mil. 2011US$)
rkna Capital stock at constant 2011 national prices (in mil. 2011US$)
rtfpna TFP at constant national prices (2011=1)
rwtfpna Welfare-relevant TFP at constant national prices (2011=1)
labsh Share of labour compensation in GDP at current national prices
delta Average depreciation rate of the capital stock
xr Exchange rate, national currency/USD (market+estimated)
pl_con Price level of CCON (PPP/XR), price level of USA GDPo in 2011=1
pl_da Price level of CDA (PPP/XR), price level of USA GDPo in 2011=1
pl_gdpo Price level of CGDPo (PPP/XR), price level of USA GDPo in 2011=1
i_cig 0/1/2, see note i_cig
i_xm 0/1/2, see note i_xm
i_xr 0/1: the exchange rate is market-based (0) or estimated (1)
i_outlier 0/1, see note i_outlier
cor_exp Correlation between expenditure shares, see note cor_exp
statcap Statistical capacity indicator (source: World Bank, developing countries only)
csh_c Share of household consumption at current PPPs
csh_i Share of gross capital formation at current PPPs
csh_g Share of government consumption at current PPPs
csh_x Share of merchandise exports at current PPPs
csh_m Share of merchandise imports at current PPPs
csh_r Share of residual trade and GDP statistical discrepancy at current PPPs
pl_c Price level of household consumption, price level of USA GDPo in 2011=1
pl_i Price level of capital formation, price level of USA GDPo in 2011=1
pl_g Price level of government consumption, price level of USA GDPo in 2011=1
pl_x Price level of exports, price level of USA GDPo in 2011=1
pl_m Price level of imports, price level of USA GDPo in 2011=1
pl_k Price level of the capital stock, price level of USA 2011=1

このセッションでは，主に以下の3つの変数を使います。

codebook[codebook$variable %in% c("country", "year", "rgdpo"), ]
## # A tibble: 3 x 2
##   variable label                                                 
##   <chr>    <chr>                                                 
## 1 country  Country name                                          
## 2 year     Year                                                  
## 3 rgdpo    Output-side real GDP at chained PPPs (in mil. 2011US$)

なお，%in% は左辺の要素が右辺のベクトルに含まれているかをチェックします。（\(a \in V\) の判定）

"A" %in% letters
## [1] FALSE
"a" %in% letters
## [1] TRUE

purrr の map_*()

パイプ

PWT のようなデータがあったとき，（順番は前後するかもしれませんが）次のような操作がよくあるパターンです。

  1. 不要な列を落とす（例えば，消費のデータは今はいらない，といった理由で）
  2. 不要な行を落とす（例えば，1960年より前のデータはいらない，といった理由で）
  3. 既存の列から計算された新しい情報を追加する
  4. グループごとに集計する
  5. 集計結果に応じて並び替える

このような処理は1つのコマンドで実行するのは難しいので，複数のステップに分けるのですが，各ステップのつながりが分かりやすいようにしておくと後々の管理が楽になります。

tidyverse ではパイプ %>% を使って複数ステップの接続を行います。

オブジェクト %>% 
  関数1() %>% 
  関数2()

という形が基本のパターンで，これは

関数2(関数1(オブジェクト))

と同じ意味です。最終出力は関数で変換されたオブジェクトです。

関数はパラメータ付きで

オブジェクト %>% 
  関数1(p = foo) %>% 
  関数2(q = bar)

のようにすることもできます。これは，

関数2(関数1(オブジェクト, p = foo), q = bar)

と同じ意味です。コードが非常に読みやすくなりますね。最初に挙げた基本ワークフローは

データ %>% 
  列を選ぶ(必要な列1, 必要な列2) %>% 
  行を取り出す(条件1, 条件2) %>% 
  列を追加する(新しい列名 = 計算式) %>% 
  グループ化する(グループ化のルール) %>% 
  集計する(集計値の名称 = 集計ルール) %>% 
  並び替える(並び替えルール)

といった形で書かれることになります。データから出発して変換操作を次々と繋いでいくイメージです。石油やガスがパイプラインを通って運ばれていく様子に似ているので，このような操作を パイプ処理 と呼びます。

tidyverse の主要な関数は

ように設計されています。パイプを使った処理がしやすいように工夫されているのです。

簡単な例

パイプ自体はmagrittr という独立したパッケーじで定義されているので，tibble でなくても使うことができます。試しに簡単な例を見てみましょう。

次のコマンドは sum(1:10) と同じ意味です。

1:10 %>% sum()
## [1] 55

1:10 の和をとってから，平方根を取る，という操作は順序どおりに関数を並べれば出来上がりです。

1:10 %>% sum() %>% sqrt()     # sqrt(sum(1:10)) と同じ 
## [1] 7.416198

最初の引数を %>% の左側のオブジェクトは関数の第1引数になります。2つ目以降の引数を指定することができます。

runif(50) %>% round(digits = 1) %>% plot(type = "l")

2つ目以降の引数に左辺のオブジェクトを渡したい場合にはドット（.）を使います。

sample(letters, 10) %>% paste0("group-", .)
##  [1] "group-h" "group-y" "group-e" "group-k" "group-q" "group-j" "group-l"
##  [8] "group-b" "group-c" "group-m"

これでパイプを使う準備ができました。

dplyr の動詞

1つの tibble を変形して何らかの結果を得るためにはデータ操作に使う関数を知っておく必要があります。中でも一番基本的なのは，select(), filter(), mutate(), group_by(), summarize(), arrange() です。

select

列を選択するのは select() です。必要な列のみを取り出すために使います。 例えば，実質GDP の時系列をプロットしたいだけなら，PWTで必要なのは

  • country
  • year
  • rgdpo

だけです。この場合は，select() を使って他の変数を落としてしまいます。

pwt90 %>% 
  select(country, year, rgdpo)
## # A tibble: 11,830 x 3
##    country  year rgdpo
##    <chr>   <dbl> <dbl>
##  1 Aruba    1950    NA
##  2 Aruba    1951    NA
##  3 Aruba    1952    NA
##  4 Aruba    1953    NA
##  5 Aruba    1954    NA
##  6 Aruba    1955    NA
##  7 Aruba    1956    NA
##  8 Aruba    1957    NA
##  9 Aruba    1958    NA
## 10 Aruba    1959    NA
## # ... with 11,820 more rows

もし，選ぶのではなく除外したいという場合には，列名にマイナス符号を付けます。

pwt90 %>% 
  select(country, year, rgdpo) %>%   # ここまでは上と同じ3列の tibble
  select(-year)                      # year を除外する
## # A tibble: 11,830 x 2
##    country rgdpo
##    <chr>   <dbl>
##  1 Aruba      NA
##  2 Aruba      NA
##  3 Aruba      NA
##  4 Aruba      NA
##  5 Aruba      NA
##  6 Aruba      NA
##  7 Aruba      NA
##  8 Aruba      NA
##  9 Aruba      NA
## 10 Aruba      NA
## # ... with 11,820 more rows

filter

列を選ぶのは select() です。では行を選ぶには？filter() を使います。 例えば，日本のデータだけを取ってきたいとしましょう。次のようにします。

pwt90 %>% 
  filter(country == "Japan")
## # A tibble: 65 x 47
##    countrycode country currency_unit  year  rgdpe  rgdpo   pop   emp   avh
##    <chr>       <chr>   <chr>         <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl>
##  1 JPN         Japan   Yen            1950 2.18e5 2.10e5  83.3  38.2 2131.
##  2 JPN         Japan   Yen            1951 2.44e5 2.33e5  84.6  39.3 2113.
##  3 JPN         Japan   Yen            1952 2.66e5 2.57e5  85.9  40.5 2094.
##  4 JPN         Japan   Yen            1953 2.79e5 2.68e5  87.1  41.8 2077.
##  5 JPN         Japan   Yen            1954 2.96e5 2.82e5  88.2  42.3 2108.
##  6 JPN         Japan   Yen            1955 3.22e5 3.05e5  89.3  43.7 2135.
##  7 JPN         Japan   Yen            1956 3.44e5 3.27e5  90.2  44.5 2188.
##  8 JPN         Japan   Yen            1957 3.68e5 3.51e5  91.0  45.7 2212.
##  9 JPN         Japan   Yen            1958 4.00e5 3.75e5  91.8  45.9 2239.
## 10 JPN         Japan   Yen            1959 4.40e5 4.11e5  92.7  46.3 2263.
## # ... with 55 more rows, and 38 more variables: hc <dbl>, ccon <dbl>,
## #   cda <dbl>, cgdpe <dbl>, cgdpo <dbl>, ck <dbl>, ctfp <dbl>,
## #   cwtfp <dbl>, rgdpna <dbl>, rconna <dbl>, rdana <dbl>, rkna <dbl>,
## #   rtfpna <dbl>, rwtfpna <dbl>, labsh <dbl>, delta <dbl>, xr <dbl>,
## #   pl_con <dbl>, pl_da <dbl>, pl_gdpo <dbl>, i_cig <dbl+lbl>,
## #   i_xm <dbl+lbl>, i_xr <dbl+lbl>, i_outlier <dbl+lbl>, cor_exp <dbl>,
## #   statcap <dbl>, csh_c <dbl>, csh_i <dbl>, csh_g <dbl>, csh_x <dbl>,
## #   csh_m <dbl>, csh_r <dbl>, pl_c <dbl>, pl_i <dbl>, pl_g <dbl>,
## #   pl_x <dbl>, pl_m <dbl>, pl_k <dbl>

日本とアメリカだけなら？プロンプトの出力では少し分かりにくいので View() と組み合わせます。アメリカのデータが残っているのが分かるでしょうか。

pwt90 %>% 
  filter(country %in% c("Japan", "United States")) %>% 
  View()

複数条件を組み合わせることもできます。AND条件はコンマで区切るだけです。(& も使えます)

pwt90 %>% 
  filter(country %in% c("Japan", "United States"), year > 2011)
## # A tibble: 6 x 47
##   countrycode country currency_unit  year  rgdpe  rgdpo   pop   emp   avh
##   <chr>       <chr>   <chr>         <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl>
## 1 JPN         Japan   Yen            2012 4.45e6 4.47e6  127.  64.3 1745 
## 2 JPN         Japan   Yen            2013 4.50e6 4.50e6  127.  64.6 1734 
## 3 JPN         Japan   Yen            2014 4.48e6 4.51e6  127.  65.0 1729 
## 4 USA         United… US Dollar      2012 1.60e7 1.59e7  315. 145.  1754.
## 5 USA         United… US Dollar      2013 1.63e7 1.62e7  317. 146.  1759.
## 6 USA         United… US Dollar      2014 1.67e7 1.66e7  319. 148.  1765.
## # ... with 38 more variables: hc <dbl>, ccon <dbl>, cda <dbl>,
## #   cgdpe <dbl>, cgdpo <dbl>, ck <dbl>, ctfp <dbl>, cwtfp <dbl>,
## #   rgdpna <dbl>, rconna <dbl>, rdana <dbl>, rkna <dbl>, rtfpna <dbl>,
## #   rwtfpna <dbl>, labsh <dbl>, delta <dbl>, xr <dbl>, pl_con <dbl>,
## #   pl_da <dbl>, pl_gdpo <dbl>, i_cig <dbl+lbl>, i_xm <dbl+lbl>,
## #   i_xr <dbl+lbl>, i_outlier <dbl+lbl>, cor_exp <dbl>, statcap <dbl>,
## #   csh_c <dbl>, csh_i <dbl>, csh_g <dbl>, csh_x <dbl>, csh_m <dbl>,
## #   csh_r <dbl>, pl_c <dbl>, pl_i <dbl>, pl_g <dbl>, pl_x <dbl>,
## #   pl_m <dbl>, pl_k <dbl>

OR 条件は | を使います。次のフィルターは「『日本かアメリカまたは人口が10億人以上』かつ『2013年以降』」のデータを抽出しています。

pwt90 %>% 
  filter(country %in% c("Japan", "United States") | pop > 1000, year > 2012)
## # A tibble: 8 x 47
##   countrycode country currency_unit  year  rgdpe  rgdpo   pop   emp   avh
##   <chr>       <chr>   <chr>         <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl>
## 1 CHN         China   Yuan Renminbi  2013 1.59e7 1.59e7 1363. 793.    NA 
## 2 CHN         China   Yuan Renminbi  2014 1.71e7 1.71e7 1369. 798.    NA 
## 3 IND         India   Indian Rupee   2013 6.36e6 6.59e6 1279. 501.  2162.
## 4 IND         India   Indian Rupee   2014 6.77e6 7.06e6 1295. 510.  2162.
## 5 JPN         Japan   Yen            2013 4.50e6 4.50e6  127.  64.6 1734 
## 6 JPN         Japan   Yen            2014 4.48e6 4.51e6  127.  65.0 1729 
## 7 USA         United… US Dollar      2013 1.63e7 1.62e7  317. 146.  1759.
## 8 USA         United… US Dollar      2014 1.67e7 1.66e7  319. 148.  1765.
## # ... with 38 more variables: hc <dbl>, ccon <dbl>, cda <dbl>,
## #   cgdpe <dbl>, cgdpo <dbl>, ck <dbl>, ctfp <dbl>, cwtfp <dbl>,
## #   rgdpna <dbl>, rconna <dbl>, rdana <dbl>, rkna <dbl>, rtfpna <dbl>,
## #   rwtfpna <dbl>, labsh <dbl>, delta <dbl>, xr <dbl>, pl_con <dbl>,
## #   pl_da <dbl>, pl_gdpo <dbl>, i_cig <dbl+lbl>, i_xm <dbl+lbl>,
## #   i_xr <dbl+lbl>, i_outlier <dbl+lbl>, cor_exp <dbl>, statcap <dbl>,
## #   csh_c <dbl>, csh_i <dbl>, csh_g <dbl>, csh_x <dbl>, csh_m <dbl>,
## #   csh_r <dbl>, pl_c <dbl>, pl_i <dbl>, pl_g <dbl>, pl_x <dbl>,
## #   pl_m <dbl>, pl_k <dbl>

最後に select() との組み合わせです。

pwt90 %>% 
  filter(year > 1960, country %in% c("Japan", "United States", "Republic of Korea")) %>% 
  select(country, year, rgdpo, pop)
## # A tibble: 162 x 4
##    country  year    rgdpo   pop
##    <chr>   <dbl>    <dbl> <dbl>
##  1 Japan    1961  517845.  94.4
##  2 Japan    1962  566475.  95.2
##  3 Japan    1963  612924.  96.2
##  4 Japan    1964  682741.  97.2
##  5 Japan    1965  724052.  98.3
##  6 Japan    1966  798912   99.2
##  7 Japan    1967  883256. 100. 
##  8 Japan    1968  995241. 101. 
##  9 Japan    1969 1117313. 103. 
## 10 Japan    1970 1233349  104. 
## # ... with 152 more rows

Non-Standard Evaluation (NSE)

pwt90 %>% 
  filter(year > 2012, country == "Japan") %>% 
  select(country, year, rgdpo, pop)
## # A tibble: 2 x 4
##   country  year    rgdpo   pop
##   <chr>   <dbl>    <dbl> <dbl>
## 1 Japan    2013 4504978.  127.
## 2 Japan    2014 4509603   127.

このコマンドでは，列名をあたかもオブジェクト名のように扱っています。year や country や rgdpo という名前はどこにも定義されていませんので，「object not found」と怒られそうなものですが，上のコマンドはうまく動きます。

実際，次のコマンドは失敗します。

year > 2012
## Error in eval(expr, envir, enclos): object 'year' not found

R の関数に渡される引数は

  • 式を評価した値だけを使う
  • 式そのものも使う

という2つのパターンがあります。「式そのもの」を使うというのは，例えば，次のコマンドの結果は，X軸，Y軸のラベルは特に指定しない限り引数として設定した式になります。もちろん，プロットする点を決めるためには「式の値」が使われています。

xval <- 1:10
random <- rnorm(10)
plot(xval, xval + random)

要するに，plot() 関数は引数を式ごと受け取って必要になったら値を調べるということをやっています。このような振る舞いは，Non-Standard Evaluation (NSE) と呼ばれ，対話的な統計環境としてのRの利便性を高めています。

dplyr の関数の多くも NSE を使っています。

  • 呼び出したときに名前が定義されていなくてもいったん式を受け取る
  • 関数が実行されているときに，tibble の中にその名前を探す
pwt90 %>% 
  select(nonexistent)
## Error in .f(.x[[i]], ...): object 'nonexistent' not found

mutate / transmute

列を新たに追加するのは mutate() を使います。

例えば，次のコマンド, 最後の行は1人あたり実質GDPを追加します。

pwt90 %>% 
  filter(year > 1960, country %in% c("Japan", "United States", "Republic of Korea")) %>% 
  select(country, year, rgdpo, pop) %>% 
  mutate(rgdpo_pc = rgdpo / pop)
## # A tibble: 162 x 5
##    country  year    rgdpo   pop rgdpo_pc
##    <chr>   <dbl>    <dbl> <dbl>    <dbl>
##  1 Japan    1961  517845.  94.4    5488.
##  2 Japan    1962  566475.  95.2    5947.
##  3 Japan    1963  612924.  96.2    6370.
##  4 Japan    1964  682741.  97.2    7022.
##  5 Japan    1965  724052.  98.3    7367.
##  6 Japan    1966  798912   99.2    8055.
##  7 Japan    1967  883256. 100.     8814.
##  8 Japan    1968  995241. 101.     9821.
##  9 Japan    1969 1117313. 103.    10896.
## 10 Japan    1970 1233349  104.    11893.
## # ... with 152 more rows

ほしいデータが rgdpo_pc だけで，rgdpo, pop が必要なくなるようなケースも考えられます。その場合，transmute() を使うことができます。これは，select() と mutate() を同時に実行してくれます。

pwt90 %>% 
  filter(year > 1960, country %in% c("Japan", "United States", "Republic of Korea")) %>% 
  transmute(country, year, rgdpo_pc = rgdpo / pop)
## # A tibble: 162 x 3
##    country  year rgdpo_pc
##    <chr>   <dbl>    <dbl>
##  1 Japan    1961    5488.
##  2 Japan    1962    5947.
##  3 Japan    1963    6370.
##  4 Japan    1964    7022.
##  5 Japan    1965    7367.
##  6 Japan    1966    8055.
##  7 Japan    1967    8814.
##  8 Japan    1968    9821.
##  9 Japan    1969   10896.
## 10 Japan    1970   11893.
## # ... with 152 more rows

脱線: ggplot2

ここまで来るとグラフを描いてみたくなります。すこし脱線して ggplot2 を紹介しましょう。

まずは，コードを見てみましょう。3行目までは前のコードと同じです。4行目からが可視化のためのコードです。

pwt90 %>% 
  filter(year > 1960, country %in% c("Japan", "United States", "Republic of Korea")) %>% 
  transmute(country, year, rgdpo_pc = rgdpo / pop) %>%    # ここまでは前と同じ
  ggplot(aes(x = year, y = rgdpo_pc, color = country)) +
    geom_line() + 
    xlab("Year") + 
    ylab("Real GDP per capita")

ggplot() という関数は「ggplot2 の可視化を始めますよ」というコマンドです。第1引数としてデータを取るので，パイプを使って上手くデータを渡すことができます。

ggplot() で描画エリアを初期化した後は，+ の記号（パイプ%>%ではなく！）を使って，グラフのレイヤー（層）を重ねていきます。aes() というのはデータを外観上の特徴に対応付けるための関数です。この例では

  • year はグラフの横軸の値を定める
  • rgdpo_pc はグラフの縦軸の値を定める
  • country はグラフの色を定める（“Japan” なら赤，など）

ということを設定しています。geom_line() というのは，折れ線グラフのレイヤーを作成するための関数です。引数が空の場合は，ggplot() に渡された データと外観特性を使って描写します。

xlab(), ylab() は軸のラベルのレイヤーを作ります。

group_by

上のような図を見ると「各国の平均成長率はどれぐらいだったのか？」といった疑問が湧いてきます。このような計算をするためには，group_by() と summarize() を使います。まずは，group_by() です。

pwt90 %>% 
  filter(year > 1960, country %in% c("Japan", "United States", "Republic of Korea")) %>% 
  transmute(country, year, rgdpo_pc = rgdpo / pop) %>% 
  group_by(country)
## # A tibble: 162 x 3
## # Groups:   country [3]
##    country  year rgdpo_pc
##    <chr>   <dbl>    <dbl>
##  1 Japan    1961    5488.
##  2 Japan    1962    5947.
##  3 Japan    1963    6370.
##  4 Japan    1964    7022.
##  5 Japan    1965    7367.
##  6 Japan    1966    8055.
##  7 Japan    1967    8814.
##  8 Japan    1968    9821.
##  9 Japan    1969   10896.
## 10 Japan    1970   11893.
## # ... with 152 more rows

見た目上なんの違いもないように見えますが，summarize() を使うと違いがはっきりと分かります。

summarize

summarize() は複数行にまたがる情報を集約して1つの結果を得るためのものです。例えば，何行あるかとか，平均はいくらかといった計算をするときに使います。通常の mean() や行数を数える n() という関数が使えます。

pwt90 %>% 
  filter(year > 1960, country %in% c("Japan", "United States", "Republic of Korea")) %>% 
  transmute(country, year, rgdpo_pc = rgdpo / pop) %>% 
  summarize(mean = mean(rgdpo_pc), n = n())
## # A tibble: 1 x 2
##     mean     n
##    <dbl> <int>
## 1 24060.   162

おっと，国ごとの結果にならなかったですね。これは group_by() を忘れたからです。

pwt90 %>% 
  filter(year > 1960, country %in% c("Japan", "United States", "Republic of Korea")) %>% 
  transmute(country, year, rgdpo_pc = rgdpo / pop) %>% 
　group_by(country) %>% 
  summarize(mean = mean(rgdpo_pc), n = n())
## # A tibble: 3 x 3
##   country             mean     n
##   <chr>              <dbl> <int>
## 1 Japan             23358.    54
## 2 Republic of Korea 13594.    54
## 3 United States     35228.    54

さて，問題の平均成長率ですが，これを計算するには

\[ Y \text{の成長率}(t\to t+1) = \log Y_{t+1} - \log Y_t \] すなわち

\[ Y \text{の平均成長率}(1 \to T) = \sum_{t=1}^{T-1} \frac{\log Y_{t+1} - \log Y_t}{T - 1} = \frac{\log Y_T - \log Y_1}{T - 1} \] という公式を使います。

pwt90 %>% 
  filter(year > 1960, country %in% c("Japan", "United States", "Republic of Korea")) %>% 
  transmute(country, year, rgdpo_pc = rgdpo / pop, ln_rgdpo_pc = log(rgdpo_pc)) %>% 
　group_by(country) %>% 
  summarize(avg_growth = mean(ln_rgdpo_pc - lag(ln_rgdpo_pc), na.rm = TRUE))
## # A tibble: 3 x 2
##   country           avg_growth
##   <chr>                  <dbl>
## 1 Japan                 0.0353
## 2 Republic of Korea     0.0633
## 3 United States         0.0203

lag() はデータを1期ずらす関数です。

arrange

上の例では行が3つだけなので大きな問題にはなりませんが，「平均成長率が最も高い国はどこか」ということを知りたい場合には，データが整列している方が何かと都合がよいです。そのために使うのが arrange() です。

pwt90 %>% 
  filter(year > 1960, country %in% c("Japan", "United States", "Republic of Korea")) %>% 
  transmute(country, year, rgdpo_pc = rgdpo / pop, ln_rgdpo_pc = log(rgdpo_pc)) %>% 
　group_by(country) %>% 
  summarize(avg_growth = mean(ln_rgdpo_pc - lag(ln_rgdpo_pc), na.rm = TRUE)) %>% # ココまでは上と同じ
  arrange(avg_growth)
## # A tibble: 3 x 2
##   country           avg_growth
##   <chr>                  <dbl>
## 1 United States         0.0203
## 2 Japan                 0.0353
## 3 Republic of Korea     0.0633

arrange() のデフォルトは小さい順です。大きい順にするには desc() を噛ませます。

pwt90 %>% 
  filter(year > 1960, country %in% c("Japan", "United States", "Republic of Korea")) %>% 
  transmute(country, year, rgdpo_pc = rgdpo / pop, ln_rgdpo_pc = log(rgdpo_pc)) %>% 
　group_by(country) %>% 
  summarize(avg_growth = mean(ln_rgdpo_pc - lag(ln_rgdpo_pc), na.rm = TRUE)) %>% # ココまでは上と同じ
  arrange(desc(avg_growth))
## # A tibble: 3 x 2
##   country           avg_growth
##   <chr>                  <dbl>
## 1 Republic of Korea     0.0633
## 2 Japan                 0.0353
## 3 United States         0.0203