Pythonで統計処理をしているが、Rでコードを書いたことがない方々のために、30分から1時間で読めるR入門を書きました。この記事の R Markdown 文書と実行環境は、こちらの GitHubレポジトリ にあります。
目次
コードブロックの左上にPythonとあるのはPythonのコード、RとあるのはRのコード、無印で各行が ## から始まるブロックは直上のブロックの実行結果の出力です。以下長文なので、ですます調ではなく、である調で書きます。
準備
本文書の GitHubレポジトリ に、インストール手順があるのでそれに従う。
準備ができたら、以下のコードを RStudio の Console に順に入力して実行する。Rスクリプトを実行するときは、拡張子 .R のファイルにコードを書いて、シェルから Rscript を起動する。
プログラミング言語Rの基本
基本型
基本型としてとりあえず、論理型 (logical : Pythonのbool)、数値 (numeric)、文字列 (character) を知っておけばよい。characterは文字列で、文字型は無い。文字コードは Unicode、エンコーディングは UTF-8 一択にする。Shift_JIS はもう忘れよう。
論理型の、真は TRUE、偽は FALSE という定数である。全部大文字である。整数と真理値を足すとき、 TRUE は1に、 FALSE は0として扱う。Python の 2 + True と、R の 2 + TRUE はどちらも 3 になる。
まずPythonの実行結果を示す。
2 + True
## 3
次にRの実行結果を示す。出力の先頭に [1] がついているが、これについては後述する。RもPythonと同様、行の # 以降は行末までコメント扱いされて実行しない。
2 + TRUE
## [1] 3
NA (欠測値 : missing value) と、変数に対応するオブジェクトが無い (NULL) ことを区別する。変数がNAかどうかは is.na で、NULLかどうかは is.null で調べる。これらの細かい仕様は省略する。
is.na(0)
## [1] FALSE
is.na(NA)
## [1] TRUE
is.null(0)
## [1] FALSE
is.null(NULL)
## [1] TRUE
一般的にNAとの計算結果はNAになる。
2 + NA
## [1] NA
2 > NA
## [1] NA
以下の結果は、C++の boost::tribool に近い。
TRUE & NA
## [1] NA
TRUE | NA
## [1] TRUE
FALSE & NA
## [1] FALSE
FALSE | NA
## [1] NA
NAは logical, numeric, character と、それぞれの型に用意されている。実は logical → numeric → character という変換は、必要なら暗黙に行われるので、あまり気にしなくてよい。さっきの 2 + TRUE がそうだった。
数値が NA (欠測値)、正の無限大 Inf、負の無限大 -Inf、NaN (非数 : Not a Number)であることを区別する。is.na, is.infinite, is.nan を使って調べる。 == では比較しないことに注意する。RはNAとNaNを区別するので便利である(入力値が欠測なことと、計算結果がNaNであることは別なので)。
is.infinite(0)
## [1] FALSE
is.infinite(Inf)
## [1] TRUE
Inf > 0
## [1] TRUE
is.infinite(-Inf)
## [1] TRUE
-Inf < 0
## [1] TRUE
is.nan(0)
## [1] FALSE
is.nan(NaN)
## [1] TRUE
関数の説明は、?関数名で得られる。RStudio の Console から ?is.na と入力すると、 is.na の説明が表示される。ネットで探すよりも速い。
?is.na
Rにスカラー変数はない。C++の std::vector てきな集合、つまり型が等しい複数の要素を、インデックスで参照できる集合が、ベクトルとして用意されている。なのでスカラー変数に見えるものは、要素が1個のベクトルである。変数への代入は <- で行う。Getsと発音する。Pythonと同様に = でも代入できるが、 <- の方が見やすいと私は思う。
n <- 2
n * 3
## [1] 6
複数の即値からベクトルを作るには、 c(15, 26, 37, 48) というように c を使う。c は combine の略である。ベクトルというか配列を作るのに [] を使うプログラミング言語が多いのでちょっと変わっているが、すぐ慣れるだろう。
score_vec <- c(15, 26, 37, 48)
name_vec <- c("foo", "bar", "poi")
score_vec
## [1] 15 26 37 48
name_vec
## [1] "foo" "bar" "poi"
出力の [1] は、この行はベクトルの1番目から順番に表示している、という意味である。さっきから [1] が必ず表示されていたのは、要素が1個のベクトルを表示していたからである。要素がたくさんあるベクトルを表示すると行を折り返すので、そのときは1より大きい数字になる。
ベクトルに限らず、配列らしくインデックスがついているデータ(PandasのDataFrameてきなものも含む)の要素は、1始まりで数える。つまりベクトル v の先頭の要素は、 v[1] で読み書きする。PythonやC++は0始まりだが、Rは異なる(FortranやStanも1始まり)。
score_vec[1]
## [1] 15
name_vec[2]
## [1] "bar"
NumPy と同様、複数のインデックスを指定することもできる。
score_vec[c(1, 4)]
## [1] 15 48
負の値を指定すると、指定したインデックス以外を選ぶ。つまり除外するインデックスを指定することができる。
score_vec[c(-2, -3)]
## [1] 15 48
0を指定すると何も選ばない。PythonやC++の癖で、先頭要素を0にするとこうなる。敢えて0を含めることで条件分岐を省けることがあるかもしれない。
score_vec[0]
## numeric(0)
ベクトル同士はcで結合する。
c(c(1, 2, 3), c(4, 5))
## [1] 1 2 3 4 5
演算
Rにスカラー変数はなく、ベクトルはあるので、演算はベクトル同士の演算になる。
足し算 c(10, 20, 30) + c(4, 5, 6) は、要素ごとの和になる。これは要素が1個のベクトル 10 + 4 = 14 の拡張と考えれば素直だろう。
c(10, 20, 30) + c(4, 5, 6)
## [1] 14 25 36
10 + 4
## [1] 14
max は、ベクトルの要素のうち、もっとも大きい値を返す。 max(c(10, 20, 30)) は c(30) 、つまり要素が1個のベクトルである。さっきの足し算が map なら、 max は reduce と考えると解りやすいだろう。
max(c(10, 20, 30))
## [1] 30
2つのベクトルの要素ごとの max を取りたい、つまり map したいなら、 pmax を使う。 pmax は長さ(要素数)nである単一または複数のベクトルの、第1要素同士を比べて最も大きい値, 第2要素同士を比べて最も大きい値 … , 第n要素同士を比べて最も大きい値を返す。
pmax(c(10, 2, 30), c(4, 50, 6))
## [1] 10 50 30
NumPy てきな broadcasting を、Rでも行うことができる。 c(10, 20, 30) + 4 は c(14, 24, 34) になる。左辺のそれぞれに4を足すというふるまいは直観的だろう。実際には右辺ベクトルの要素を、左辺ベクトルの要素数と同じになるまで繰り返している (repeatという)。
c(10, 20, 30) + 4
## [1] 14 24 34
ちなみに同じ値を繰り返したければ、 rep を使う。
rep(4, 3)
## [1] 4 4 4
Assertion
Pythonでは、 () をつけずに assert を書く。
assert 2 ** 4 == 16
assert 2 ** 4 != 16
## Error in py_call_impl(callable, dots$args, dots$keywords): AssertionError:
Rでは assertthat を使う。パッケージを使う前に import するのは、Pythonと同様である。
library(assertthat)
assertthat::assert_that((2**4) == 16)
## [1] TRUE
assertthat::assert_that((2**4) != 16)
## Error: (2^4) not not equal to 16
ベクトルの要素数とrange
要素数は NROW で調べる。 length もあるが忘れていい。
NROW(c(10, 20, 30))
## [1] 3
NROW(4)
## [1] 1
NROW(c())
## [1] 0
Python の range てきな連番を作ろう。要素数n個のベクトルvに対して 1, 2, …, n という連番の整数は、 seq_len(n) で作る。
seq_len(3)
## [1] 1 2 3
1:3 とも書ける。
1:3
## [1] 1 2 3
seq_len(NROW()) でベクトルのインデックスを作れる。seq_len(NROW(c(10, 20, 30))) は要素が3個なので c(1, 2, 3) になる。 seq_len(NROW(c())) は空ベクトルが返る。
seq_len(NROW(c(10, 20, 30)))
## [1] 1 2 3
seq_len(NROW(c()))
## integer(0)
続きは等差数列の項で説明する。
リスト
型が混在する集合は list に格納する。Pythonの list やC言語の構造体と似ている。以下のコードに “” (二重引用符) が無い箇所があることに注意する。
a_lst <- list(name = "foo", score = 80, year = 2019)
a_lst
## $name
## [1] "foo"
##
## $score
## [1] 80
##
## $year
## [1] 2019
NROW(a_lst)
## [1] 3
names(a_lst)
## [1] "name" "score" "year"
要素には $属性名 、 [[属性名]] 、インデックスでアクセスできる。Pythonのenumerableな感じだ。
a_lst$score
## [1] 80
a_lst[["score"]]
## [1] 80
a_lst[[2]]
## [1] 80
一重カッコ [属性名] は、要素ではなくリストを返す。 is.list はオブジェクトがリストかどうか返す。
a_lst["score"]
## $score
## [1] 80
a_lst[2]
## $score
## [1] 80
is.list(a_lst["score"])
## [1] TRUE
is.list(a_lst[["score"]])
## [1] FALSE
ベクトルと同様に、インデックスを指定して部分リストを得ることができる。属性名を指定してもよい。
a_lst[c(1, 2)]
## $name
## [1] "foo"
##
## $score
## [1] 80
a_lst[c("name", "score")]
## $name
## [1] "foo"
##
## $score
## [1] 80
リスト同士はcで結合する。Pythonもそうだが、リストを長くするのか、リストにリストを入れて入れ子にするのかで、コードを使い分ける必要がある。
こちらはリストを長くする。
list_1 <- list("foo")
list_2 <- list(80)
list_3 <- list(2019)
list_3_elements <- c(list_1, list_2, list_3)
list_3_elements
## [[1]]
## [1] "foo"
##
## [[2]]
## [1] 80
##
## [[3]]
## [1] 2019
NROW(list_3_elements)
## [1] 3
こちらはリストを入れ子にする。
nested_list <- list(first = list_1, second = list_2)
nested_list
## $first
## $first[[1]]
## [1] "foo"
##
##
## $second
## $second[[1]]
## [1] 80
要素の型がすべて同じリストは、 unlist でベクトルに変換することができる。逆も然りで、 as.list を使ってベクトルをリストにする。後で説明するが、Rでは一般的に as.型名 で型変換を行う。
identical は、二つのオブジェクトが同値ならTRUE、違っていたらFALSEを返す。identicalといいつつ、同一オブジェクト(雑に言うと同じメモリ番地に置かれている)かどうかではなく、値が同じかどうかを調べる。
original_vec <- 1:5
list_from_vec <- as.list(original_vec)
converted_vec <- unlist(list_from_vec)
assertthat::assert_that(identical(original_vec, converted_vec))
## [1] TRUE
DataFrame (tibble)
library(tidyverse)
Pandas の DataFrame と同様に、Rにも data.frame がある (Rが先だったかもしれないが)。だが tibble の方が便利なのでそちらを使う。まず説明用の架空データを作る。
df <- tibble(name = c("foo", "bar"), score = c(80, 60), year = c(2019, 2020))
| name | score | year |
|---|---|---|
| foo | 80 | 2019 |
| bar | 60 | 2020 |
変数を print すると今作った tibble を文字列として表示する。というより実は、変数名だけ入力すると print して表示していたのだった。Python の REPL と同じである。
print(df)
## # A tibble: 2 x 3
## name score year
## <chr> <dbl> <dbl>
## 1 foo 80 2019
## 2 bar 60 2020
View するとExcelっぽいウィンドウに表示する。
View(df)
tibble は、縦に標本を並べ、横に標本の属性を並べる。ざっくりいうと、縦に長いベクトル=すべての標本についてのそれぞれの属性を、横にいろいろと並べたものがtibbleと思えばよい。なので、先のリストの動作がそのまま当てはまる。
まず score のベクトルを取得する。Pandas がそうであるように、tibble は連想配列らしく使える。ベクトルが返るので Pandas の Series とは少し異なる。
df$score
## [1] 80 60
df[["score"]]
## [1] 80 60
df[[2]]
## [1] 80 60
score列だけを切り出した tibble を取得する。
df["score"]
## # A tibble: 2 x 1
## score
## <dbl>
## 1 80
## 2 60
df[2]
## # A tibble: 2 x 1
## score
## <dbl>
## 1 80
## 2 60
行=標本も指定しよう。この記法は NumPy っぽく、Pandas のように loc とは書かない。
df[[1, "score"]]
## [1] 80
df[1, "score"]
## # A tibble: 1 x 1
## score
## <dbl>
## 1 80
行=標本を指定して全属性を取り出すときは、 , の後を空白にする。Pandasのように : は書かない。
df[1, ]
## # A tibble: 1 x 3
## name score year
## <chr> <dbl> <dbl>
## 1 foo 80 2019
ここでは一列ずつ、一行ずつ取り出したが、リストと同様に複数の要素を取り出すこともできる。
Matrix
Matrix は tibble に比べれば覚えることが少ないので、最低限のことだけ説明する。
tibble と異なり、すべての要素を同じ型にする。まずベクトルから行列を作る。 byrow で並べ方を指定する。
matrix(data = 1:6, nrow = 2, byrow = FALSE)
## [,1] [,2] [,3]
## [1,] 1 3 5
## [2,] 2 4 6
matrix(data = 1:6, nrow = 2, byrow = TRUE)
## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] 4 5 6
要素がすべて型の同じ tibble から行列を作る。
df_num <- tibble(p = 1:2, q = 3:4, r = 5:6)
mat_num <- as.matrix(df_num)
mat_num
## p q r
## [1,] 1 3 5
## [2,] 2 4 6
as_tibbleでtibbleに変換する。
df_mat_num <- as_tibble(mat_num)
df_mat_num
## # A tibble: 2 x 3
## p q r
## <int> <int> <int>
## 1 1 3 5
## 2 2 4 6
assertthat::assert_that(identical(df_num, df_mat_num))
## [1] TRUE
演算、転置、3次元以上の配列 (array) の生成などはもちろんできるので、必要になったら調べればよい。
関数
Rの関数は、Pythonのラムダ式である。
Pythonでは以下のように、引数の2乗を返す無名関数を定義して、それを my_square という変数に格納し、こう呼び出せる。
my_square = lambda n: n * n
my_square(9)
## 81
Rでは function と書く。
my_square <- function(n) {
n * n
}
my_square(9)
## [1] 81
Pythonと同様に、引数にキーワードをつけてもよい。引数が多いときは、キーワードをつけるとソースコードが読みやすく保守しやすくなるだろう。
my_square(n = 9)
## [1] 81
Pythonと同様に、引数のキーワードありとなしを混在することもできる。まずキーワードがある引数を当てはめ、残りのパラメータをまだ値が無い先頭のパラメータから順に当てはめる。
my_pow <- function(base, index) {
base**index
}
my_pow(index = 2, 3)
## [1] 9
do.call を使って、Pythonのunpack, Rubyのsplatのように、listを引数に展開することもできる。
args <- list(base = 3, index = 2)
do.call(my_pow, args)
## [1] 9
さて、Rの関数はすべて式(expression)である。関数が最後に評価したこと(雑に言うと {} 内の最後の行)が関数の返り値になる。関数定義の途中で処理を終えるのでなければ、 return と明示する必要ない(この辺はRubyっぽい)。
Rで負の数の対数を取ると、警告を出してNaNを返すが、代わりに-Infを返す関数を作ろう。ifも式なのでこう書ける。
my_log10 <- function(x) {
if (x <= 0) {
-Inf
} else {
log10(x)
}
}
log10(-1)
## Warning: 計算結果が NaN になりました
## [1] NaN
my_log10(-1)
## [1] -Inf
my_log10(0)
## [1] -Inf
my_log10(0.1)
## [1] -1
my_log10(10)
## [1] 1
最初に述べた通り、Rにスカラー変数はないので、この引数xはベクトルである。ならば、ベクトルを受け取ってベクトルを返すmap動作がRらしいだろう。Pythonの Pythonic に対応する言葉が、R-like である。上記のコードがどうなるか試してみよう。
my_log10(c(-1, 0, 0.1, 10))
## Warning in if (x <= 0) {: 条件が長さが 2 以上なので、最初の 1 つだけが使われます
## [1] -Inf
上手くいってないようだ。改善方法は二つある。一つ目の方法は、 if ではなく ifelse を使って、ベクトルの要素ごとにthenな値とelseな値を選ぶことだ。ただしこれは、log10(負の値)が警告を出すので、見た目がよろしくない。警告ではなくエラーになる関数では使えないだろう。
my_log10_alt1 <- function(xs) {
ifelse(xs <= 0, -Inf, log10(xs))
}
my_log10_alt1(c(-1, 0, 0.1, 10))
## Warning in ifelse(xs <= 0, -Inf, log10(xs)): 計算結果が NaN になりました
## [1] -Inf -Inf -1 1
二つ目の方法は、Pythonの list comprehension てきにmapすることだ。このようにmapを使いこなすとよいだろう。 purrr::map_dbl は要素が数値のベクトルを返す。
my_log10_alt2 <- function(xs) {
purrr::map_dbl(xs, function(x) {
if (x <= 0) {
-Inf
} else {
log10(x)
}
})
}
my_log10_alt2(c(-1, 0, 0.1, 10))
## [1] -Inf -Inf -1 1
Rでfor文は使わないので、忘れてよい。Pythonでfor文を使う処理があったら、常にmapとreduceを使うことを検討しよう。mapの例の次に、reduceの使い方としてフィボナッチ数列を示そう。0と負の数についてはエラーになるが、それでいいだろう。
my_fib <- function(xs) {
purrr::map_dbl(xs, function(x) {
tail(purrr::reduce(.x = seq_len(x), .init = integer(), function(acc, i) {
if (i == 1) {
c(1)
} else if (i == 2) {
c(1, 1)
} else {
c(acc[2], sum(acc))
}
}), 1)
})
}
my_fib(1:10)
## [1] 1 1 2 3 5 8 13 21 34 55
Pythonで外の lexical scope にある変数を変更するには、nonlocal を使う。以下のコードでは、outer() 呼び出しは1を返す。
def outer():
z = 0
def inner():
nonlocal z
z = z + 1
return -2
inner()
return z
outer()
## 1
Rでは <<- を使って、 lexical scope の外にある変数に代入する。
outer <- function() {
z <- 0
inner <- function() {
z <<- z + 1
-2
}
inner()
z
}
outer()
## [1] 1
PythonもRも、lexical scopeの範囲は、細かい説明を省いてざっくりいうと関数(defやfunction)である。なので関数が入れ子になっている場合は、nonlocalが必要になることがある。C++は if (cond) {} や何もキーワードがつかない {} も自動変数の範囲と寿命を決めているが、Rはそうではないので注意する(つまり if などは lexical scope を変えない)。
いろいろな処理と込み入った話題
文字列処理
Rでは文字列の整形をするには、昔ながらの sprintf を使う。
sprintf("地点=%s, 温度=%.1f", "横浜", 11.8)
## [1] "地点=横浜, 温度=11.8"
文字列を連結するときは、 paste を使う。 sep 引数が文字列の区切り文字で、デフォルトは半角空白1個である。英語で英単語を連結して文を作るときはこれでよいが、日本語では困ることがあるだろう。区切り文字を空文字列にするには sep 引数で明示的に指定するか、 paste0 を使う。Rの文字列は + 演算子で連結できない。
paste("神奈川県", "横浜市", "中区")
## [1] "神奈川県 横浜市 中区"
paste("神奈川県", "横浜市", "中区", sep = " ")
## [1] "神奈川県 横浜市 中区"
paste0("神奈川県", "横浜市", "中区")
## [1] "神奈川県横浜市中区"
paste("神奈川県", "横浜市", "中区", sep = ":")
## [1] "神奈川県:横浜市:中区"
以下のコードは意外にも、要素が3個のベクトルを返す。 paste はベクトルをベクトルに map する関数である。
words <- c("神奈川県", "横浜市", "中区")
paste(words)
## [1] "神奈川県" "横浜市" "中区"
collapse を指定すると、文字列を区切り文字で連結 (join) して、単一の文字列を返す reduce な動作になる。
paste(words, sep = "", collapse = ":")
## [1] "神奈川県:横浜市:中区"
より高度な文字列処理は stringrパッケージを使う。PythonもRもたくさんの関数があるが、代表的なものとして正規表現で文字列を分解する例を示そう。Python の reモジュールは full matching 、R の stringrパッケージは partial matching である。
import re
re.match(r".{2}(\D+)(\d+)", "01Year2022").groups()
## ('Year', '2022')
stringr::str_match("01Year2022", "(\\D+)(\\d+)")
## [,1] [,2] [,3]
## [1,] "Year2022" "Year" "2022"
強制型変換
既に書いたが、logical, numeric, character の右側への変換は、必要なら暗黙に行われる。 paste 関数で試そう。書式を指定しなければ、よきに計らった文字列表記になる。
paste("Text", exp(1), 2, FALSE)
## [1] "Text 2.71828182845905 2 FALSE"
リストとベクタの変換で既に示したように、型変換は一般に as.型名 でできる。Rは変数名に . を含めることができ、PythonやC++でインスタンスメソッドの前に . をつけるときのような特別な意味は . にはない。
as.numeric("-2.7")
## [1] -2.7
as.integer("-2.7")
## [1] -2
as.character(-2.7)
## [1] "-2.7"
無限大と文字列に書くことができる。Pythonの float() 呼び出しと同じである。こちらがPythonで
float("Inf")
## inf
float("-infinity")
## -inf
こちらがRである。
as.numeric("Inf")
## [1] Inf
as.numeric("-infinity")
## [1] -Inf
型変換できないときは、警告を出してNAが返る。Pythonと違って、エラーにはならない。
float("abc")
## Error in py_call_impl(callable, dots$args, dots$keywords): ValueError: could not convert string to float: 'abc'
as.numeric("abc")
## Warning: 強制変換により NA が生成されました
## [1] NA
NAは後でどうにかするので警告を無視したいなら、 suppressWarnings を使う。
suppressWarnings(as.numeric("abc"))
## [1] NA
等差数列
Pythonの range と同様のことをしたければ、Rには seq とか 1:n とかがある。Pythonと異なりRは閉区間であり、区間の開始と終了の両方を要素に含む。また、開始番号が終了番号より小さいと、1ずつ減る数列が得られる。ベクトルのインデックスを得る場合に、これだと空のベクトルのインデックスが c(1, 0) になってまずい。 seq_len を積極的に使おう。
seq(4)
## [1] 1 2 3 4
1:4
## [1] 1 2 3 4
seq_len(4)
## [1] 1 2 3 4
seq(0)
## [1] 1 0
1:0
## [1] 1 0
seq_len(0)
## integer(0)
seq(-2)
## [1] 1 0 -1 -2
1:-2
## [1] 1 0 -1 -2
seq_len(-2)
## Error in seq_len(-2): 引数は非負の整数に変換できなければなりません
さて差が整数ではない等差数列は、丸め誤差が発生するので、予想した長さにならないことがある。以下はPythonのコードである。 numpy.linspace は数列の長さを指定するので正しい個数になるが、 numpy.arrange は丸め誤差の影響で個数が合わないことがある。
import numpy as np
for x in range(50, 53):
print("{} {}".format(len(np.linspace(0.0, 1.0, x)),
len(np.arange(0.0, 1.0, 1.0/(x - 1)))))
## 50 50
## 51 50
## 52 51
同様のことを、Rでは seq の length.out 引数と、 by 引数で区別する。この例ではPythonと違って丸め誤差の影響を受けなかったが、両引数を使い分ける必要がある。
purrr::map_chr(50:52, function(x) {
sprintf(
"%d %d", NROW(seq(from = 0.0, to = 1.0, length.out = x)),
NROW(seq(from = 0.0, to = 1.0, by = 1.0 / (x - 1)))
)
})
## [1] "50 50" "51 51" "52 52"
これまで明言しなかったが、関数呼び出しの結果や、変数を評価した結果はコンソールに表示される。そのままだと大量の出力でコンソールが埋まって困ることがあるだろう。関数の結果を変数で受け止めるのが一つの方法だが、もう一つの方法は関数呼び出しの結果を invisible で受け止めることである。
いかにもわざとらしい例だが、以下の例は invisible があると何も表示しないし、 invisible がなければ最大20000個(通常は1000個で打ち切り)を表示しようとする。
invisible(1:20000)
集合演算
二つのベクトルA, Bについて、和集合(A, Bの少なくとも片方に属する要素)、積集合(A, B の両方に属する要素)、差集合(A=左にあって、B=右にない要素)を、それぞれ union, intersect, setdiff で求める。 setdiff は非対称なので、引数を入れ替えると結果が変わることがある。
union(1:4, 3:7)
## [1] 1 2 3 4 5 6 7
intersect(1:4, 3:7)
## [1] 3 4
setdiff(1:4, 3:7)
## [1] 1 2
setdiff(3:7, 1:4)
## [1] 5 6 7
ベクトルAの各要素がベクトルBに含まれているかどうかを調べるために、 A %in% B を使うと、含まれている=TRUE または 含まれていない=FALSE からなるベクトルを返す。演算子のヘルプを見るには、 ?`%in%` というように演算子を ` (backtick) で囲む。
1:7 %in% 3:4
## [1] FALSE FALSE TRUE TRUE FALSE FALSE FALSE
3:4 %in% 1:7
## [1] TRUE TRUE
日時と時刻
日時(date)と時刻(timestamp)を上手く扱うのは難しいので、よくできたパッケージを上手に使いこなす。一年の初めからどれだけの割合(0以上1以下)が経過したか指定すると、その日時を表示する関数を作る。うるう年は考慮するがうるう秒は考慮しない。あるツイートにあった、1月4日のお昼には、もう2021年の「1%」が経過しているのは本当だろうか?
library(lubridate)
elapsed_time_in_year <- function(year, ratio) {
epoch <- paste0(year, "-01-01 00:00:00")
x <- lubridate::ymd_hms(epoch, tz = "Asia/Tokyo")
x + lubridate::dseconds((365 + lubridate::leap_year(x)) * 24 * 60 * 60 * ratio)
}
elapsed_time_in_year(year = 2021, ratio = 0.01)
## [1] "2021-01-04 15:36:00 JST"
elapsed_time_in_year(year = 2020, ratio = 0.01)
## [1] "2020-01-04 15:50:24 JST"
elapsed_time_in_year(year = 2021, ratio = 0.162)
## [1] "2021-03-01 03:07:12 JST"
elapsed_time_in_year(year = 2020, ratio = 0.162)
## [1] "2020-02-29 07:00:28 JST"
時間 (duration)を説明するついでに関数合成しよう。 functional::Compose は二つの関数を合成する。ここでは hh:mm:ss 書式の時間の文字列を受け取って、何秒か返す関数を作る。 function を使っても同じことをできるが、引数を書いて値を転送しなくても済むのですっきりする。
library(functional)
hms_to_seconds <- functional::Compose(lubridate::hms, lubridate::seconds)
hms_to_seconds("00:01:04")
## [1] "64S"
hms_to_seconds("01:00:04")
## [1] "3604S"
クラス
PythonやC++てきなクラス、つまりクラスにメンバ関数とメンバ変数を定義したデータ構造をRでも作れる。まずPythonでカウンターを作る例を示す。
class Counter:
"""A classic example of classes"""
__count = 0
def __init__(self):
"""Initialize member variable(s)"""
self.__count = 0
def increment(self):
"""Increment the count"""
self.__count = self.__count + 1
def get(self):
"""Get the count"""
return self.__count
ct = Counter()
print(ct.get())
## 0
ct.increment()
ct.increment()
print(ct.get())
## 2
Rでは R6パッケージを使う。ほぼPythonと同様に書ける。public, private と明記するところがC++みたいだ。
library(R6)
Counter <- R6::R6Class("Counter",
public = list(
# Initialize member variable(s)
initialize = function() {
private$count <- 0
},
# Increment the count
increment = function() {
private$count <- private$count + 1
},
# Get the count
get = function() {
private$count
}
),
private = list(
count = 0
)
)
ct <- Counter$new()
print(ct$get())
## [1] 0
ct$increment()
ct$increment()
print(ct$get())
## [1] 2
参照渡しと copy-on-modify
Rでは、呼び出された関数は引数を変更できるが、その変更は呼び出し側には伝わらない(そうでない場合についてはすぐ後で述べる)。これはPythonとの大きな違いである。より正確に言うと、引数を変更すると copy-on-modify される。引数を変更しなければコピーするコストが発生しないので、実質的には低コストな参照渡しのように動作する。
Pythonで、配列の要素をすべて2倍にするとき、こう書ける。
import numpy as np
def doubler(x):
for i, item in enumerate(x):
x[i] = item * 2
print(np.sum(x))
assert np.sum(x) == 20
v = np.ones(10)
doubler(v)
## 20.0
v
## array([2., 2., 2., 2., 2., 2., 2., 2., 2., 2.])
Rで同様のコードを書くとどうなるだろうか? 忘れていたfor文を敢えて使って書こう。
doubler <- function(x) {
for (i in seq_len(NROW(x))) {
x[i] <- x[i] * 2
}
print(sum(x))
assertthat::assert_that(sum(x) == 20)
}
v <- rep(1, 10)
doubler(v)
## [1] 20
## [1] TRUE
v
## [1] 1 1 1 1 1 1 1 1 1 1
doubler の中では二倍になったようだが(和が20なので)、呼び出し元の変数vは変わらなかった。これが copy-on-modify である。
Rでは引数を上書きすることは諦めて、値を返す関数を素直にmapで書こう。 ~ .x * 2 は無名関数の略記法である。今は深入りしないが、いずれ役に立つことがあるだろう。
doubler <- function(x) {
purrr::map_dbl(x, ~ .x * 2)
}
doubler(rep(1, 10))
## [1] 2 2 2 2 2 2 2 2 2 2
クラスのインスタンス(R6オブジェクト)は参照渡しである。呼び出された関数が引数に変更を加えると、引数を渡した呼び出し側にも変更が反映される。
increment_in_callee <- function(obj) {
obj$increment()
print(obj$get())
}
ct_caller <- Counter$new()
invisible(increment_in_callee(ct_caller))
## [1] 1
print(ct_caller$get())
## [1] 1
Pythonも参照渡しなので、参照渡しの方がなじみがあるだろう。
def increment_in_callee(obj):
obj.increment()
print(obj.get())
ct_caller = Counter()
increment_in_callee(ct_caller)
## 1
print(ct_caller.get())
## 1
デフォルト値
Pythonで関数呼び出しに失敗した場合にデフォルト値を得ようとしたら、try-except構文を使うだろう。
import math
def str_to_float(x: str) -> float:
value = math.nan
try:
value = float(x)
except ValueError:
pass
return(value)
math.isnan(str_to_float("1.25"))
## False
math.isnan(str_to_float("x.y.z"))
## True
Rでは purrr::possibly を使って、デフォルト値を返す関数を作る。 stop はエラーを発生させる関数である。
func_may_unsupported <- function(x) {
stop("Unsupported")
}
safe_func <- purrr::possibly(func_may_unsupported, NA)
safe_func(1)
## [1] NA
変数が存在するかどうかは、Pythonでは、 locals と globals で調べると分かる。
a_glogal_var = 1
def check_vars_exists():
a_local_var = 2
assert "a_glogal_var" in globals()
assert "a_local_var" in locals()
assert "missing_var" not in locals()
check_vars_exists()
Rでは exists で分かる。
a_glogal_var <- 1
check_vars_exists <- function() {
a_local_var <- 2
assertthat::assert_that(exists("a_glogal_var")) &
assertthat::assert_that(exists("a_local_var")) &
assertthat::assert_that(!exists("missing_var"))
}
check_vars_exists()
## [1] TRUE
Pythonと同様にデフォルト引数を書くこともできる。だがデフォルト引数がいつ評価されるかは、PythonでもRでも厄介な問題を生むので、個人的にはお勧めしない。以下はPythonでありがちなコードである。
def add(x=[]):
x.append(1)
return x
list_a = add([])
print(list_a)
## [1]
add(list_a)
## [1, 1]
list_b = add([])
print(list_b)
## [1]
add(list_b)
## [1, 1]
おそらく以下の挙動は望んでいないだろう。
add()
## [1]
add()
## [1, 1]
それはさておき、Rのデフォルト引数も、見た目はPythonと同様である。
real_to_complex <- function(real, imaginary = 0.0) {
list(r = real, im = imaginary)
}
real_to_complex(real = 1)
## $r
## [1] 1
##
## $im
## [1] 0
real_to_complex(real = 1, imaginary = 2)
## $r
## [1] 1
##
## $im
## [1] 2
引数が無いことは missing で分かる。NULLを指定することと、引数が無いことは異なる。この例ならデフォルト引数を使っても大差ないが、引数が無いことが分かるとできることがあるかもしれない。
real_to_complex_alt <- function(real, imaginary) {
im_part <- if (missing(imaginary)) {
0
} else {
imaginary
}
list(r = real, im = im_part)
}
real_to_complex_alt(real = 1)
## $r
## [1] 1
##
## $im
## [1] 0
real_to_complex_alt(real = 1, imaginary = 2)
## $r
## [1] 1
##
## $im
## [1] 2
real_to_complex_alt(real = 1, imaginary = NULL)
## $r
## [1] 1
##
## $im
## NULL
関数のデフォルト引数の値を調べるために、Pythonでは inspect.signature を使う。Rでは formals を使う。以下は正規表現で文字列を分割する関数のデフォルト引数である。
import inspect
inspect.signature(re.split)
## <Signature (pattern, string, maxsplit=0, flags=0)>
formals(stringr::str_split)
## $string
##
##
## $pattern
##
##
## $n
## [1] Inf
##
## $simplify
## [1] FALSE
Assertionで実行を止める
Runtime assertion で条件を満たさないときに、処理を止めることができる。浮動小数が丸め誤差を除いて大体同じかどうかを、 assertthat::are_equal で許容誤差 tol (toleranceの略) を与えて比べる。
assertthat::are_equal(x = 6.0, y = 6.1, tol = 0.5)
## [1] TRUE
assertthat::are_equal は二つの数が大体同じかどうかを返すが、同じという条件を満たさなくても処理は止まらない。単にFALSEを返すだけである。
check_if_equal <- function(value1, value2) {
assertthat::are_equal(x = value1, y = value2, tol = 0.5)
}
check_if_equal(value1 = -5, value2 = 6)
## [1] FALSE
条件を満たさないときに処理を止めたければ、 assertthat::assert_that を併用する。
stop_if_not_equal <- function(value1, value2) {
assertthat::assert_that(assertthat::are_equal(x = value1, y = value2, tol = 0.5))
}
stop_if_not_equal(value1 = -5, value2 = 6)
## Error: value1 not equal to value2
コマンドライン引数を解析する
コマンドライン引数は、 commandArgs で文字列ベクトルとして取得する。 trailingOnly=FALSE にするとすべての引数、 trailingOnly=TRUE にすると –args および以降の引数(–argsがなければ空文字列ベクトル)を取得する。R処理系に渡す引数を無視して、Rスクリプトだけに渡す引数を取得するために、 trailingOnly=TRUE を使うとよい。
RStudio 上でRスクリプトを実行したときは、コマンドライン引数としてオプションを渡すことができない。対策としては、グローバル変数(仮に g_args と名付ける)があったらコマンドライン引数の代わりに扱い、そうでなければ実際のコマンドライン引数を commandArgs を呼び出して取得する、とすればよい。変数があるかどうかは、先ほどの exists を使う。
get_commandline_args <- function() {
if (exists("g_args")) {
g_args
} else {
commandArgs(trailingOnly = TRUE)
}
}
get_commandline_args()
## character(0)
g_args <- c("param1", "param2")
get_commandline_args()
## [1] "param1" "param2"
CSVファイルを読んで集計する
Rの基本が分かったので、CSVファイルをRで読み込んで描画しよう。Rにはサンプルデータがあってそれを使うのが定番であるが、そういう描画例はウェブで多数見つかるのと、Pythonで同じデータを描画できることを目的に、横浜市のオープンデータを使用する(詳細は データ出典 )。
ここではダウンロードしたデータが、 incoming_yokohama/ サブディレクトリに置いてあるものとする。ウェブ上のデータは更新されているので、実行時期によって結果が変わるかもしれない。
| 内容とリンク先 | 更新日 | ファイル名 |
|---|---|---|
| 人口と世帯数の推移 | 最新版(2020年まで) | jinkosetai-sui.csv |
| 男女別人口及び世帯数－行政区 | 令和3年12月(2021年12月) | e1yokohama2112.csv |
ディレクトリを作る
CSVファイルを読んで集計した結果を出力するディレクトリを作る。既にディレクトリがあるならそのままにする。
out_dirname <- "output"
dir.create(out_dirname, showWarnings = FALSE)
assertthat::assert_that(dir.exists(out_dirname))
## [1] TRUE
CSVファイルを読む
readr::read_csv を使う。
library(tidyverse)
df_timeseries <- readr::read_csv("incoming_yokohama/jinkosetai-sui.csv")
## Rows: 132 Columns: 4
## -- Column specification -------------------------------------------------------------------------
## Delimiter: ","
## chr (1): 年（和暦）
## dbl (3): 年（西暦）, 世帯数［世帯］, 人口［人］
##
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.
| 年（和暦） | 年（西暦） | 世帯数［世帯］ | 人口［人］ |
|---|---|---|---|
| 明治22 | 1889 | 27209 | 121985 |
| 明治23 | 1890 | 27835 | 127987 |
| 明治24 | 1891 | 29070 | 132627 |
| 明治25 | 1892 | 29269 | 143252 |
| 明治26 | 1893 | 29942 | 152142 |
| 明治27 | 1894 | 29974 | 160439 |
さて print で表示すると、列名に `` がついているのが分かる。
print(df_timeseries)
## # A tibble: 132 x 4
## `年（和暦）` `年（西暦）` `世帯数［世帯］` `人口［人］`
## <chr> <dbl> <dbl> <dbl>
## 1 明治22 1889 27209 121985
## 2 明治23 1890 27835 127987
## 3 明治24 1891 29070 132627
## 4 明治25 1892 29269 143252
## 5 明治26 1893 29942 152142
## 6 明治27 1894 29974 160439
## 7 明治28 1895 30124 170252
## 8 明治29 1896 30474 179502
## 9 明治30 1897 31584 187453
## 10 明治31 1898 31765 191251
## # ... with 122 more rows
これは列名に ( を含むので、 $ の後に列名を書けないということである。
df_timeseries$年（和暦）
## Error: <text>:1:17: 想定外の入力です
## 1: df_timeseries$年（
## ^
という風に、 ( が上手く構文解析できないのでエラーになってしまった。構文中に埋め込む名前として使えるものを syntactic name 、そうでないものを non-syntactic name と呼ぶ。詳しい解説は Advanced R にあるが、ここでは間にあわせの回避策だけ述べよう。
一つ目は、列名を ` (backtick)で囲むことである。
df_timeseries$`年（和暦）`
## [1] "明治22" "明治23" "明治24" "明治25" "明治26" "明治27" "明治28" "明治29" "明治30" "明治31"
## [11] "明治32" "明治33" "明治34" "明治35" "明治36" "明治37" "明治38" "明治39" "明治40" "明治41"
## [21] "明治42" "明治43" "明治44" "大正1" "大正2" "大正3" "大正4" "大正5" "大正6" "大正7"
## [31] "大正8" "大正9" "大正10" "大正11" "大正12" "大正13" "大正14" "大正15" "昭和2" "昭和3"
## [41] "昭和4" "昭和5" "昭和6" "昭和7" "昭和8" "昭和9" "昭和10" "昭和11" "昭和12" "昭和13"
## [51] "昭和14" "昭和15" "昭和16" "昭和17" "昭和18" "昭和19" "昭和20" "昭和21" "昭和22" "昭和23"
## [61] "昭和24" "昭和25" "昭和26" "昭和27" "昭和28" "昭和29" "昭和30" "昭和31" "昭和32" "昭和33"
## [71] "昭和34" "昭和35" "昭和36" "昭和37" "昭和38" "昭和39" "昭和40" "昭和41" "昭和42" "昭和43"
## [81] "昭和44" "昭和45" "昭和46" "昭和47" "昭和48" "昭和49" "昭和50" "昭和51" "昭和52" "昭和53"
## [91] "昭和54" "昭和55" "昭和56" "昭和57" "昭和58" "昭和59" "昭和60" "昭和61" "昭和62" "昭和63"
## [101] "平成1" "平成2" "平成3" "平成4" "平成5" "平成6" "平成7" "平成8" "平成9" "平成10"
## [111] "平成11" "平成12" "平成13" "平成14" "平成15" "平成16" "平成17" "平成18" "平成19" "平成20"
## [121] "平成21" "平成22" "平成23" "平成24" "平成25" "平成26" "平成27" "平成28" "平成29" "平成30"
## [131] "令和1" "令和2"
二つ目は、列名を colnames で付け替えることである。もちろん入力するデータを確認して、どんな形式か分かった上で、それぞれの列に対応する名前をつける。
colnames(df_timeseries) <- c("JPyear", "Year", "Household", "Population")
| JPyear | Year | Household | Population |
|---|---|---|---|
| 明治22 | 1889 | 27209 | 121985 |
| 明治23 | 1890 | 27835 | 127987 |
| 明治24 | 1891 | 29070 | 132627 |
| 明治25 | 1892 | 29269 | 143252 |
| 明治26 | 1893 | 29942 | 152142 |
| 明治27 | 1894 | 29974 | 160439 |
グラフを描く
まず単純なグラフを描こう。横軸を年(西暦)、縦軸を人口の折れ線グラフを描く。たった3行で済む。
g <- ggplot(df_timeseries)
g <- g + geom_line(aes(x = Year, y = Population))
plot(g)
いろいろと見た目を変えたり、要素を付け加えることもできる。詳しくは R for Data Scienceに例がたくさんある。
- X軸の刻みを10年間隔にする
- フォントを Segoe UI にして(フォントがなければデフォルトのsansフォントにする)、フォントのサイズを指定する
- 折れ線の色を青くして太くする
- タイトルをつける
- 背景を白にする
- 高さを幅の0.6倍にする
year_unit <- 10
x_min <- ceiling(min(df_timeseries$Year) / year_unit) * year_unit
x_max <- floor(max(df_timeseries$Year) / year_unit) * year_unit
x_breaks <- unique(c(
min(df_timeseries$Year),
seq(from = x_min, to = x_max, by = year_unit), max(df_timeseries$Year)
))
x_head <- head(x_breaks, 2)
x_tail <- tail(x_breaks, 2)
x_breaks <- setdiff(
x_breaks,
c(x_head[2 * (diff(x_head) < year_unit)], x_tail[1 * (diff(x_tail) < year_unit)])
)
library(extrafont)
font_name <- "Segoe UI"
font_set <- extrafont::fonts()
if (is.null(font_set) || !(font_name %in% font_set)) {
font_name <- "sans"
}
g <- ggplot(df_timeseries)
g <- g + geom_line(aes(x = Year, y = Population), color = "navy", size = 2)
g <- g + scale_x_continuous(breaks = x_breaks)
g <- g + ggtitle("Yokohama's Population")
g <- g + theme_bw()
g <- g + theme(
aspect.ratio = 0.6,
text = element_text(family = 16),
axis.text = element_text(family = font_name, size = 12),
axis.title = element_text(family = font_name, size = 16),
plot.title = element_text(family = font_name, size = 20)
)
plot(g)
tibbleの列名は英語にしても、図の表示は日本語にしたいということがあるだろう。 Migu 1M フォントで表示する。
font_name_jp <- "Migu 1M"
g <- ggplot(df_timeseries)
g <- g + geom_line(aes(x = Year, y = Population), color = "navy", size = 2)
g <- g + scale_x_continuous(breaks = x_breaks)
g <- g + ggtitle("横浜市の人口")
g <- g + xlab("年(西暦)")
g <- g + ylab("人口")
g <- g + theme_bw()
g <- g + theme(
aspect.ratio = 0.6,
text = element_text(family = font_name_jp),
axis.text = element_text(family = font_name_jp, size = 12),
axis.title = element_text(family = font_name_jp, size = 16),
plot.title = element_text(family = font_name_jp, size = 20)
)
plot(g)
plotly::ggplotly を使うと、インタラクティブに拡大縮小したり、線にカーソルを置いて各年の人口を表示できたりする。ただし、出力されるHTMLファイルが数Mbytes増えるのと、PDF化できなくなる。コメントアウトを外して試してほしい。
library(plotly)
## plotly::ggplotly(g)
図をPNGファイルに保存するには、 ggsave を使う。 file.path はディレクトリ名 out_dirname と ファイル名を結合する。Pythonで os.path.join を使うのと同様に、 / で結合するよりも正しい作法である。
png_filename <- file.path(out_dirname, "yokohama_population.png")
ggsave(png_filename, plot = g, dpi = 160)
図をPNGファイルに保存するために png と dev.off で囲んでもよいが、その場合は高さと幅のピクセル数を固定する必要がある。ピクセル数を指定できるのは便利だが、縦横比を変えたくないときは ggsave の方がよいだろう。
png(filename = png_filename, width = 800, height = 600)
g <- ggplot(df_timeseries)
g <- g + geom_line(aes(x = Year, y = Population))
plot(g)
dev.off()
## png
## 2
CSVファイルを加工する
Pandasてきに、CSVファイルを加工する方法を説明する。横浜市各区の人口を読み込む。その後、列名を英語にする。
df_population_data <- readr::read_csv("incoming_yokohama/e1yokohama2112.csv")
## Rows: 19 Columns: 12
## -- Column specification -------------------------------------------------------------------------
## Delimiter: ","
## chr (1): 市区名
## dbl (10): 全国地方公共団体コード, 面積［平方キロメートル］, 世帯数［世帯］, 人口総数［人］, 男［人］, 女［人］, １世帯当たり人員［人］, 人口密度［人／平方キロメ...
## date (1): 年月日
##
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.
print(colnames(df_population_data))
## [1] "年月日" "全国地方公共団体コード"
## [3] "市区名" "面積［平方キロメートル］"
## [5] "世帯数［世帯］" "人口総数［人］"
## [7] "男［人］" "女［人］"
## [9] "１世帯当たり人員［人］" "人口密度［人／平方キロメートル］"
## [11] "届出による前月比増減の世帯数［世帯］" "届出による前月比増減の人口［人］"
colnames(df_population_data) <- c(
"date", "code", "ward", "area",
"household", "population", "male", "female",
"per_household", "density",
"diff_household", "diff_population"
)
今後使いそうな列だけ残して、後は削除する。列を選ぶには dplyr::select を使う。 %>% は、コマンドラインというかシェルのパイプと同様、前段で処理したデータを後段に流す。この書き方に慣れるとRが好きになるだろう。実際には前段の返り値を後段の第一引数にしており、Pandasのメソッドチェーン(method chaining)と同じである。
df_population <- df_population_data %>%
dplyr::select(c("ward", "area", "household", "population", "per_household", "density"))
| ward | area | household | population | per_household | density |
|---|---|---|---|---|---|
| 横浜市 | 437.78 | 1767422 | 3773982 | 2.14 | 8621 |
| 鶴見区 | 33.23 | 145761 | 295940 | 2.03 | 8906 |
| 神奈川区 | 23.72 | 130479 | 247572 | 1.90 | 10437 |
| 西区 | 7.03 | 57297 | 104604 | 1.83 | 14880 |
| 中区 | 21.50 | 85101 | 150567 | 1.77 | 7003 |
| 南区 | 12.65 | 104535 | 197596 | 1.89 | 15620 |
| 港南区 | 19.90 | 96569 | 215244 | 2.23 | 10816 |
| 保土ケ谷区 | 21.93 | 99744 | 206963 | 2.07 | 9437 |
| 旭区 | 32.73 | 107257 | 243403 | 2.27 | 7437 |
| 磯子区 | 19.05 | 79178 | 166373 | 2.10 | 8733 |
| 金沢区 | 30.96 | 90459 | 197796 | 2.19 | 6389 |
| 港北区 | 31.40 | 177090 | 359753 | 2.03 | 11457 |
| 緑区 | 25.51 | 80277 | 183249 | 2.28 | 7183 |
| 青葉区 | 35.22 | 134217 | 311146 | 2.32 | 8834 |
| 都筑区 | 27.87 | 86286 | 214880 | 2.49 | 7710 |
| 戸塚区 | 35.79 | 123411 | 284210 | 2.30 | 7941 |
| 栄区 | 18.52 | 53331 | 120523 | 2.26 | 6508 |
| 泉区 | 23.58 | 63477 | 152107 | 2.40 | 6451 |
| 瀬谷区 | 17.17 | 52953 | 122056 | 2.30 | 7109 |
横浜市の合計と各区を分離しよう。 dplyr::filter を使って、～区とそれ以外を抜き出す。 str_end は文字列がある文字列で終わっているかどうかをTRUEまたはFALSEで返す。Pythonの endswith と同じである。
df_wards <- df_population %>%
dplyr::filter(stringr::str_ends(ward, "区"))
df_whole_city <- df_population %>%
dplyr::filter(!stringr::str_ends(ward, "区"))
もしくはこう書ける。
df_population %>% dplyr::filter(ward != "横浜市")
| ward | area | household | population | per_household | density |
|---|---|---|---|---|---|
| 横浜市 | 437.78 | 1767422 | 3773982 | 2.14 | 8621 |
dplyr::select は列名を文字列で指定したが、 dplyr::filter の ward は “” が無いので文字列ではない。これは先ほど説明したシンボル(symbol)である。今回のように各区はwardという列にある、と分かっているときはこのようにハードコーディングしてよいが、実行時に文字列変数で与えたいことがあるだろう。そのときは !!rlang::sym を使う。
library(rlang)
ward_column_name <- "ward"
df_wards_alt <- df_population %>%
dplyr::filter(!!rlang::sym(ward_column_name) != "横浜市")
assertthat::assert_that(identical(df_wards, df_wards_alt))
## [1] TRUE
rlang::expr を使って、 !!rlang::sym がどう展開されたかを確認すると、先ほどのコードと同じことが分かる。
rlang::expr(df_population %>%
dplyr::filter(!!rlang::sym(ward_column_name) != "横浜市"))
## df_population %>% dplyr::filter(ward != "横浜市")
実は dplyr 1.0.0 から across を使えるので、ここで rlang::sym を使う必要はない。後で rlang::sym を使う場面があるので、敢えてここで使った。
df_wards_alt_across <- df_population %>%
dplyr::filter(across(all_of(ward_column_name)) != "横浜市")
assertthat::assert_that(identical(df_wards_alt, df_wards_alt_across))
## [1] TRUE
実は dplyr::select は、文字列だけでなくシンボルを受け付ける。なので “” が無くてもよかった。わざとらしい例だが、複数の文字列をシンボルにするには、 rlang::syms を使う。それと、 !! ではなく !!! を使う。 !! と !!! が否定の否定、否定の否定の否定、ではなく演算子というのが珍しい。
column_names <- c("ward", "area", "household", "population", "per_household", "density")
df_population_alt <- df_population_data %>%
dplyr::select(!!!rlang::syms(column_names))
assertthat::assert_that(identical(df_population, df_population_alt))
## [1] TRUE
それぞれの行を集計する
df_wards には、 1世帯当たり人員[人] (per_household 列) と 人口密度[人/平方キロメートル] (density 列) がある。これらは、人口 (population 列)、世帯数 (household 列)、面積 (area 列) から計算できる。 dplyr::mutate で実際に計算して列を作ろう。
df_wards_calc <- df_wards %>%
dplyr::mutate(
per_household_calc = population / household,
density_calc = population / area
)
元データの per_household と計算で求めた per_household_calc 、元データの density と計算で求めた density_calc が等しいか確認しよう。浮動小数ベクトルが丸め誤差を除いて大体同じかどうかは、 assertthat::are_equal で許容誤差 tol を与えて比べる。
assertthat::assert_that(assertthat::are_equal(
x = df_wards_calc$per_household,
y = df_wards_calc$per_household_calc, tol = 0.01
))
## [1] TRUE
assertthat::assert_that(assertthat::are_equal(
x = df_wards_calc$density,
y = df_wards_calc$density_calc, tol = 1.0
))
## [1] TRUE
計算した1世帯当たり人員と人口密度があっていたようだ。
assertで条件を満たさないときに実行を停止するのではなく単に同じかどうか知りたければ、 dplyr::near で許容誤差 tol を与えて比べる。 dplyr::near は二つの要素を先頭から順に比較して、大体同じ=TRUE、違う=FALSEというベクトルを返す。全部TRUEかどうかは all で分かる。
are_near_values <- dplyr::near(
x = df_wards_calc$per_household,
y = df_wards_calc$per_household_calc, tol = 0.01
)
are_near_values
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
all(are_near_values)
## [1] TRUE
ちなみに dplyr::mutate の代入先もシンボルであり、実行時に文字列を与えるときは rlang::sym を使う。このときは = ではなく := を使う。代入先ではなく代入元がシンボルの場合も同様である。
dst_column_name <- "per_household_calc"
df_wards_calc_dst <- df_wards %>%
dplyr::mutate(!!rlang::sym(dst_column_name) := population / household)
assertthat::assert_that(
assertthat::are_equal(
df_wards_calc$per_household_calc, df_wards_calc_dst$per_household_calc
)
)
## [1] TRUE
src1_column_name <- "population"
src2_column_name <- "household"
df_wards_calc_src <- df_wards %>%
dplyr::mutate(per_household_calc :=
!!rlang::sym(src1_column_name) / !!rlang::sym(src2_column_name))
assertthat::assert_that(
assertthat::are_equal(df_wards_calc$per_household_calc, df_wards_calc_src$per_household_calc)
)
## [1] TRUE
それぞれの列を集計する
全区の世帯数と人口を足して、横浜市全体 df_whole_city と同じかどうか確かめる。 dplyr::summarize_at と sum を使って、列の全要素を足す。
df_wards_sum <- df_wards_calc %>%
dplyr::summarize_at(c("household", "population"), sum)
| household | population |
|---|---|
| 1767422 | 3773982 |
assertthat::assert_that(df_whole_city$household == df_wards_sum$household)
## [1] TRUE
assertthat::assert_that(df_whole_city$population == df_wards_sum$population)
## [1] TRUE
各区の合計は市全体と一致したようだ。
ここでは世帯数と人口を足したが、それは区の名前、1世帯当たり人員、人口密度を足しても意味が無いからである。そして面積は足していなかった。列名を挙げるのがめんどくさければ、列が数字だったら足す、という操作をしたい。そういうときは、 dplyr::summarize と across を使う。1世帯当たり人員と人口密度を除いてから、数値つまり区の名前以外の列を where(is.numeric) で選んで足す。
df_wards_sum_all <- df_wards %>%
dplyr::select(-c("per_household", "density")) %>%
dplyr::summarize(across(where(is.numeric), sum))
面積もあっているようだ。
assertthat::assert_that(assertthat::are_equal(x = df_whole_city$area, y = df_wards_sum_all$area, tol = 1.0))
## [1] TRUE
行をグループ化する
グループ化の例として、横浜市各区を 税務署の管轄地域 でまとめよう。まず税務署と区の対応表を作る。
tax_offices <- c(
"神奈川", "神奈川", "鶴見", "戸塚", "戸塚", "戸塚", "保土ケ谷",
"保土ケ谷", "保土ケ谷", "緑", "緑", "緑", "横浜中", "横浜中",
"横浜南", "横浜南", "横浜南", "横浜南"
)
ward_names <- c(
"神奈川区", "港北区", "鶴見区", "戸塚区", "栄区", "泉区",
"保土ケ谷区", "旭区", "瀬谷区", "緑区", "青葉区", "都筑区", "中区", "西区",
"南区", "磯子区", "金沢区", "港南区"
)
df_tax_offices <- tibble(tax_office = tax_offices, ward = ward_names)
区の名前で結合する。Rには連想配列がないが、tibbleで連想配列てきなテーブルを作ればよい。連想配列の要素を個々に引くより、 inner join や left join する方が手軽である。
df_ward_tax_offices <- dplyr::inner_join(df_tax_offices, df_wards, by = "ward")
税務署ごとに世帯数と人口をまとめよう。 dplyr::group_by を使って、行をグループ化できる。それぞれのグループについて、列ごとにすべての行を足す。
df_group_population <- df_ward_tax_offices %>%
dplyr::select(c("household", "population", "tax_office")) %>%
dplyr::group_by(tax_office) %>%
dplyr::summarize(across(everything(), sum)) %>%
dplyr::ungroup()
| tax_office | household | population |
|---|---|---|
| 横浜中 | 142398 | 255171 |
| 横浜南 | 370741 | 777009 |
| 戸塚 | 240219 | 556840 |
| 神奈川 | 307569 | 607325 |
| 鶴見 | 145761 | 295940 |
| 保土ケ谷 | 259954 | 572422 |
| 緑 | 300780 | 709275 |
鶴見が少ないのは鶴見区しかないからだが、横浜中も少ない。おそらく人口の割に事業所が多いからだと予想するが、実際そうなっているか確認するのは読者への課題とする。ちなみに衆議院選挙小選挙区は区や市をまたぐので、今回の目的には使えなかった。
カテゴリ変数
グループに用いた税務署名は文字列である。これをカテゴリ値にすることができる。Rではカテゴリをfactorと呼ぶ。 as.factor で文字列をカテゴリに変換する。
df_factor <- df_ward_tax_offices %>%
dplyr::mutate(tax_office = as.factor(tax_office))
カテゴリを表示すると、見た目は文字列のままだが、実際は文字列に対応する整数を振ったものだと分かる。
df_factor$tax_office
## [1] 神奈川 神奈川 鶴見 戸塚 戸塚 戸塚 保土ケ谷 保土ケ谷 保土ケ谷 緑
## [11] 緑 緑 横浜中 横浜中 横浜南 横浜南 横浜南 横浜南
## Levels: 横浜中 横浜南 戸塚 神奈川 鶴見 保土ケ谷 緑
as.integer(df_factor$tax_office)
## [1] 4 4 5 3 3 3 6 6 6 7 7 7 1 1 2 2 2 2
カテゴリの整数の番号を固定したいことがある。例えばグラフの凡例を描くとき、凡例はカテゴリの整数が小さい順(1, 2, …)に表示される。 forcats::fct_relevel を使うと、整数の順番を固定できる。
カテゴリの順番を五十音順にしよう。つまり、神奈川(かながわ)、鶴見(つるみ)、戸塚(とつか)、保土ケ谷(ほどがや)、緑(みどり)、横浜中(よこやまなか)、横浜南(よこはまみなみ)の順にする。
df_factor <- df_ward_tax_offices %>%
dplyr::mutate(tax_office = as.factor(tax_office)) %>%
dplyr::mutate(tax_office = forcats::fct_relevel(
tax_office,
c("神奈川", "鶴見", "戸塚", "保土ケ谷", "緑", "横浜中", "横浜南")
))
df_factor$tax_office
## [1] 神奈川 神奈川 鶴見 戸塚 戸塚 戸塚 保土ケ谷 保土ケ谷 保土ケ谷 緑
## [11] 緑 緑 横浜中 横浜中 横浜南 横浜南 横浜南 横浜南
## Levels: 神奈川 鶴見 戸塚 保土ケ谷 緑 横浜中 横浜南
as.integer(df_factor$tax_office)
## [1] 1 1 2 3 3 3 4 4 4 5 5 5 6 6 7 7 7 7
先ほどとは整数の割り当てが異なる。図にしたとき、凡例が五十音順になることを確かめよう。
library(RColorBrewer)
g <- ggplot(df_factor)
g <- g + geom_point(aes(x = per_household, y = density, color = tax_office), size = 4, shape = 4)
g <- g + geom_text(aes(x = per_household, y = density, label = ward),
nudge_x = -0.03, size = 2
)
g <- g + scale_color_manual(values = brewer.pal(NROW(unique(df_factor$tax_office)), "Dark2"))
g <- g + ylim(0, 16000)
g <- g + xlab("1世帯当たり人員[人]")
g <- g + ylab("人口密度[人/平方キロメートル]")
g <- g + guides(color = guide_legend(title = "税務署"))
g <- g + theme_bw()
g <- g + theme(
aspect.ratio = 0.8,
text = element_text(family = font_name_jp),
legend.position = "right",
legend.text = element_text(family = font_name_jp, size = 16),
legend.title = element_text(family = font_name_jp, size = 16),
axis.text = element_text(family = font_name_jp, size = 16),
axis.title = element_text(family = font_name_jp, size = 14),
strip.text = element_text(family = font_name_jp, size = 20),
plot.title = element_text(family = font_name_jp, size = 20)
)
plot(g)
facet_wrap(~ tax_office, nrow=2) と一行追加するだけで、カテゴリ別の表を作る。ggplot2すごい。
g <- ggplot(df_factor)
g <- g + geom_point(aes(x = per_household, y = density, color = tax_office), size = 4, shape = 4)
g <- g + geom_text(aes(x = per_household, y = density, label = ward),
nudge_x = -0.03, size = 2
)
g <- g + scale_color_manual(values = brewer.pal(NROW(unique(df_factor$tax_office)), "Dark2"))
g <- g + ylim(0, 16000)
g <- g + xlab("1世帯当たり人員[人]")
g <- g + ylab("人口密度[人/平方キロメートル]")
g <- g + guides(color = guide_legend(title = "税務署"))
g <- g + theme_bw()
g <- g + facet_wrap(~tax_office, nrow = 2)
g <- g + theme(
aspect.ratio = 0.8,
text = element_text(family = font_name_jp),
legend.position = "none",
axis.text = element_text(family = font_name_jp, size = 9),
axis.title = element_text(family = font_name_jp, size = 10),
strip.text = element_text(family = font_name_jp, size = 12),
plot.title = element_text(family = font_name_jp, size = 12)
)
plot(g)
CSVファイルを出力する
readr::write_excel_csv を使って、BOMつきUTF-8でCSVファイルを書く。こうしておくとExcelで読める (BOMがないとExcelで読めない)。Pandasと同様、gzip形式で読み書きできる。
readr::write_excel_csv(
df_ward_tax_offices,
file.path(out_dirname, "yokohama_tax_office.csv")
)
readr::write_excel_csv(
df_ward_tax_offices,
file.path(out_dirname, "yokohama_tax_office.csv.gz")
)
JSONファイルを読み書きする
CSVファイルには、それがどんなデータかというメタデータを記述しづらい。メタデータをJSONファイルに書くとよいだろう。例えば横浜市と茅ヶ崎市について、このように書く。city:漢字名、kana:かな表記、roman:ローマ字表記、designated:政令指定都市かどうか、n_wards:政令指定都市なら区の数、tax_offices:所管税務署、date:最終更新日、とする。
{
"city" : "横浜市",
"kana" : "よこはま",
"roman" : "Yokohama",
"designated" : true,
"n_wards" : 18,
"tax_offices" : ["神奈川", "鶴見", "戸塚", "保土ケ谷", "緑", "横浜中", "横浜南"],
"date": "2022/01/01"
}
{
"city" : "茅ヶ崎市",
"kana" : "ちがさき",
"roman" : "Chigasaki",
"designated" : false,
"tax_offices" : "藤沢",
"date": "2022/01/01"
}
library(jsonlite)
data_yokohama <- jsonlite::fromJSON("incoming_metadata/yokohama.json")
print(data_yokohama)
## $city
## [1] "横浜市"
##
## $kana
## [1] "よこはま"
##
## $roman
## [1] "Yokohama"
##
## $designated
## [1] TRUE
##
## $n_wards
## [1] 18
##
## $tax_offices
## [1] "神奈川" "鶴見" "戸塚" "保土ケ谷" "緑" "横浜中" "横浜南"
##
## $date
## [1] "2022/01/01"
data_chigasaki <- jsonlite::fromJSON("incoming_metadata/chigasaki.json")
print(data_chigasaki)
## $city
## [1] "茅ヶ崎市"
##
## $kana
## [1] "ちがさき"
##
## $roman
## [1] "Chigasaki"
##
## $designated
## [1] FALSE
##
## $tax_offices
## [1] "藤沢"
##
## $date
## [1] "2022/01/01"
class で変数の型が分かる。
class(data_chigasaki)
## [1] "list"
JSONファイルをリストとして読み込んだことが分かる。Rにはスカラー変数が無いので、要素 (木構造の末端)はすべてベクトルである。論理型、数値、文字列型が上手く推測されたことが分かる。日付は必要なら変換しよう。
data_yokohama$date <- lubridate::ymd(data_yokohama$date)
data_yokohama$date
## [1] "2022-01-01"
class(data_yokohama$date)
## [1] "Date"
存在しない要素を取得するとNULLが返る。エラーにはならない。茅ヶ崎市は政令指定都市ではないので、区の数はNULLである。
data_yokohama$n_wards
## [1] 18
data_chigasaki$n_wards
## NULL
JSONファイルを書き出してみる。
out_json_filename <- file.path(out_dirname, "output_yokohama.json")
jsonlite::write_json(data_yokohama, path = out_json_filename, pretty = TRUE)
出力結果はこうなる。要素が1個でも配列になるのが、いかにもRらしい。
{
"city": ["横浜市"],
"kana": ["よこはま"],
"roman": ["Yokohama"],
"designated": [true],
"n_wards": [18],
"tax_offices": ["神奈川", "鶴見", "戸塚", "保土ケ谷", "緑", "横浜中", "横浜南"],
"date": ["2022-01-01"]
}
ツール
Jupyter Notebook と R Markdown
Pythonでコードを試しながら書く場合、Jupyter Notebook を使っている方が多いだろう。Jupyter Notebook とRという組み合わせも可能である。もう一つの選択肢は、 R Markdown である。まさにこの文書の元は R Markdown で、それを Markdown に変換したもののHTML表示を見ているはずである。
コードの整形と静的解析
Pythonのコードは、 autopep8 で整形して、pylint, pep8, flake8, mypy でコードの問題点を静的解析して指摘するだろう(他にもツールがあるかもしれない)。
Rは styler パッケージでコードを整形して、 lintr パッケージでコードの問題点を静的解析して指摘するとよい。 styler はRスクリプト(.Rファイル)だけでなく、R Markdown (.Rmdファイル)も整形できる。
lintr で問題点を挙げるとき、スクリプトが完成してから lintr を通すと大量の問題点を指摘されて降参しかねないので、毎日細かく通すとよい。CI (Continuous Integration) を使っているなら、毎登録 (git push) ごとにすればよい。もちろん lintr の前に、 styler で自動的に整形することもできる。
lintr よくある指摘は以下の通りである。
- Cyclomatic complexity が大きすぎる。関数が複雑すぎて、理解もテストもできないだろう。小さな関数に分割して、Cyclomatic complexity を減らす。 pylint だと、関数が長い、分岐が多いといった類似の指摘事項が出る。
-
1:NROW(aTibble)は、aTibble が空行のときに誤動作する。代わりにseq_len(NROW(aTibble))を使う。このようによくある間違ったイディオムを指摘する。Pythonでは、range(len(anObj))の代わりにenumrate(anObj)を使え、と指摘されることがある。 - ローカル変数が使われていない。似たような変数と取り違えているか、本当に使われなくなった変数が残っているかどちらかなので対処する。変数は要らなくても、その変数を返す関数呼び出しは必要かもしれない。
- 書式の良し悪し、例えば
=の前後に空白をいれるかどうかは、 styler に任せて手作業で直さない。 - 変数名に対する指摘は、仕方ない場合もある。指摘に対応しないこともしばしばある。
コードの論理行数
cloc パッケージで、ファイルの論理行、コメント行、空行を数えることができる。Rだけでなくさまざまなプログラミング言語に対応している。Perl が必要なので、Windows の場合は perl.exe がある場所にPATHを通しておく。ファイル名だけでなくURLを指定することもできる。A Shiny app developed with the golem package にあるファイルの行数を測ろう(以下の行数は測定当時)。
library(cloc)
cloc::cloc_pkg() %>%
dplyr::select(all_of(c("language", "file_count", "loc", "blank_lines", "comment_lines"))) %>%
dplyr::filter(language %in% c("R", "C++")) %>%
to_embedded_table()
| language | file_count | loc | blank_lines | comment_lines |
|---|---|---|---|---|
| R | 41 | 1154 | 201 | 425 |
| C++ | 2 | 51 | 10 | 19 |
パッケージ化
Pythonで wheelパッケージを作成して成果物を配布できるのと同様、Rもパッケージを作成して成果物を配布できる。詳細は後述の参考文献を読んで頂くとして、RStudio の Buildタブを開いて Check ボタンを押すだけで、自作したRパッケージを検査して問題点を挙げることができる。配布するパッケージにドキュメントなどがそろっていて十分な品質かどうかを、簡単に確認できる。
次に進むには
最初にデータサイエンス100本ノック（構造化データ加工編）を解いて、Pythonの DataFrame と、Rの tibble の違いを覚えるとよいでしょう。
参考書としては、以下の順番が一案です。
- 最初に読むのは、「Rではじめるデータサイエンス」 (原著: R for Data Science )
- Rをもっと使いこなすなら、「Rクックブック 第2版」 (原著: R Cookbook 2nd Edition )
- RStudioとR Markdownの使い方を学ぶなら、「再現可能性のすゝめ ―RStudioによるデータ解析とレポート作成― (Wonderful R 3)」 または、 「改訂2版 Rユーザのための RStudio［実践］入門 ―tidyverseによるモダンな分析フローの世界―」
- 「みんなのR 第2版」 (原著: R for Everyone 2nd edition)、「パーフェクトR」もよいでしょう。
- RStudioだけでなくコマンドラインを使いこなしたいなら、「バイオインフォマティクスデータスキル」 (原著: Bioinformatics Data Skills)
- プログラミング言語Rを深く知りたいなら、「R言語徹底解説」(ただし和訳は版が古い、原著: Advanced R )。 公式マニュアルより読みやすいと思います。
- Rパッケージの作成を通じて、配布可能な成果物の作り方を学びたいなら、「Rパッケージ開発入門 テスト、文書化、コード共有の手法を学ぶ」(原著 : R Packages )
- ブラウザから対話的に表示を変えたいなら、Shinyの使い方を学ぶ。「RとShinyで作るWebアプリケーション」を読んでから、Mastering Shinyを読み、そのあとEngineering Production-Grade Shiny Appsを読むとよいでしょう。
- Rだけでは足りなくて、C++で処理を書きたいなら、 みんなのRcpp を読んで Rcppパッケージを使います。
参考文献
| 書名 | 著者 | 年 | 出版社またはリンク |
|---|---|---|---|
| R for Data Science | Garrett Grolemund and Hadley Wickham | 2017 | https://r4ds.had.co.nz/ |
| R Cookbook 2nd Edition | James (JD) Long and Paul Teetor | 2019 | https://rc2e.com/ |
| R for Everyone 2nd edition | Jared P. Lander | 2017 | Addison-Wesley Professional |
| Bioinformatics Data Skills | Vince Buffalo | 2015 | O’Reilly Media |
| Advanced R 2nd edition | Hadley Wickham | 2019 | https://adv-r.hadley.nz/ |
| R Packages | Hadley Wickham | 2015 | http://r-pkgs.had.co.nz/ |
| Mastering Shiny | Hadley Wickham | 2020 | https://mastering-shiny.org/ |
| Engineering Production-Grade Shiny Apps | Colin Fay, Sébastien Rochette, Vincent Guyader and Cervan Girard | 2021 | https://engineering-shiny.org/ |
| 書名 | 著者 | 年 | 出版社またはリンク |
|---|---|---|---|
| Rではじめるデータサイエンス | Hadley Wickham, Garrett Grolemund 著/黒川利明 訳/大橋真也 監修 | 2017 | オライリー・ジャパン |
| Rクックブック 第2版 | J.D. Long, Paul Teetor 著/大橋真也 監訳/木下哲也 訳 | 2020 | オライリー・ジャパン |
| 再現可能性のすゝめ ―RStudioによるデータ解析とレポート作成― (Wonderful R 3) | 高橋康介 著/石田基広 監修/市川太祐, 高橋康介, 高柳慎一, 福島真太朗, 松浦健太郎 編 | 2018 | 共立出版 |
| 改訂2版 Rユーザのための RStudio［実践］入門 ―tidyverseによるモダンな分析フローの世界― | 松村優哉, 湯谷啓明, 紀ノ定保礼, 前田和寛 著 | 2021 | 技術評論社 |
| みんなのR 第2版 | Jared P. Lander 著/高柳慎一, 津田真樹, 牧山幸史, 松村杏子, 簑田高志 監修 | 2018 | マイナビ出版 |
| パーフェクトR | Rサポーターズ 著 | 2017 | 技術評論社 |
| バイオインフォマティクスデータスキル オープンソースツールを使ったロバストで再現性のある研究 | Vince Buffalo 著/片山俊明, 川島秀一, 鈴木治夫, 山本泰智 監訳/酒匂寛, 山村吉信 訳 | 2020 | オライリー・ジャパン |
| R言語徹底解説 | Hadley Wickham 著/石田基広, 市川太祐, 高柳慎一, 福島真太朗 訳 | 2016 | 共立出版 |
| Rパッケージ開発入門 テスト、文書化、コード共有の手法を学ぶ | Hadley Wickham 著/瀬戸山雅人, 石井弓美子, 古畠敦 訳 | 2016 | オライリー・ジャパン |
| RとShinyで作るWebアプリケーション | 梅津雄一, 中野貴広 著 | 2018 | C&R研究所 |
| みんなのRcpp | Masaki E. Tsuda 著 | 2020 | https://teuder.github.io/rcpp4everyone_ja/ |
データ出典
本文書には、 横浜市のオープンデータ を編集・加工したものを埋め込んでいます。これらは 利用条件等 に記載の、 CC BY 4.0 に基づいています。