Rではじめるデータサイエンス | 書評

Rではじめるデータサイエンスを選んだ理由

AIの分野で用いられている2大言語は、「R」と「Python」です。
Pythonは以前触ってみたので、今回はRに触ってみました。

使った本は、「Rではじめるデータサイエンス」です。本屋で目立つところにあったので、人気があると思って選びました。

一言感想

Rの感想を一言で述べると、「奥が深すぎて、習得するのに他の言語の数倍の時間がかかりそう」です。
本自体の感想は、「データの可視化から入るなど、興味を持って進められる構成。ただし、データの型や構造の説明が後半になってしまっていて、しかも解り難いので、僕としては中途半端な感覚で結末を迎えた」です。

本文

Rは、短い関数を書くだけで、複雑な図を描画したり、複雑な統計計算をしてくれる強力なパッケージ群を備えるソフトです。自分で関数をコツコツと作るより、専門家が開発してくれた既存の関数を選んで並べるといったイメージです。
csvの読み込みや日付の変換の関数、それにパッケージの数が腐るほどあるので、覚えていればすぐに結果が得られる反面、知らなければ調べるだけで、時間がどんどん過ぎて行ってしまいます。
つまり、便利だから片手間に使おうというお気軽ソフトではなくて、使う以上は腰を据えて取り組むべきソフトだと思いました。

ところどころに練習問題があることは、理解を深めるために有益でした。
ただ、問題以前に記載されていた範囲の内容では、解けないことが多かったので、調査の時間をかなり要しました。この記事の下方に解答例を載せておきましたので、よろしければご参照ください。

また、この本を読み進める上で、コードをほぼ全て打ち込みました。
たまに記載通りにタイプしても動作しない場合がありましたので、以下に「誤植等」としてまとめました。よろしければこちらもお使いください。

「R」と「Python」に触れてみた現在、自分でAIを組むとしたらPythonを選びます。ただ、現時点ではAIを組む必要性を感じていません。

誤植等

P.65
練習問題の2.
‘ 以下1行に誤植あり(cancel ledと、lとlの間にスペースがある)
誤:not_cancel led %>% count(tailnum, wt = distance)
正:not_cancelled %>% count(tailnum, wt = distance)

P.98
図の画面を出すには、「Tools」→「Grobal Options」と辿る

P.127
誤:x = col_double(),
正:x = col_integer(),

P.183
誤:vowels = str_count(word, “[aeiou]”),
consonants = str_count(word, “[^aeiou]”)
正:vowels = str_count(words, “[aeiou]”),
consonants = str_count(words, “[^aeiou]”)

P.188
練習問題
誤:replace_all()を使って、str_to_lower()の・・・
正:str_replace_all()を使って、str_to_lower()の・・・

P.219
13.4.2
誤:・・・日や月といった「人間の」時間間隔に合わせるのです。・・・
正:・・・日や月といった「人間の」時間感覚に合わせるのです。・・・

P.294
頁の先頭に次の一行が抜けている
mu <- list(5, 10, -3)

P.339
誤:data_grid(cut, .model = mod_diamond2) %>%
正:data_grid(cut, lcarat = -0.515, color = “G”, clarity = “SI1”, .model = mod_diamond2) %>%

P.345
日本語環境の場合
誤:filter(wday == “Sat“) %>%
正:filter(wday == ““) %>%
P.346も同様に”Sat”は”土”である

P.366
誤:nest(year:gdpPercap[TK17])
正:nest(year:gdpPercap)

P.389
“fuel-economy.Rmd”は、full pathでないと読み込めなかった

P.394
誤:”Fuel efficiency generally decrease with”
“engine size”,
)
正:”Fuel efficiency generally decrease with”,
“engine size”
) ,

P.399
誤:”Increasing engine size is \nrelated to”
正:”Increasing engine size is \nrelated to”,

練習問題の解答例

練習問題の解答が無くて結構困ったので、解けたところを載せておきます。間違えていたらごめんなさい。分からない問題については、正直に分からないとしてあります。なお、この解答を使うことによって、万が一事故や損害が発生したとしても当方が責任を負うことはできません。あくまでも自己責任でご利用ください。

また、練習問題と同じ章内で、練習問題より前に書かれているコードが実行されていないと動かない場合がありますので、ご注意ください。

1章 ggplot2によるデータ可視化

1.2.2 ggplotを作る

# 1.
# 四角いだけで何も表示されない

# 2.
# 32行×11列

# 3.
# f: 前輪駆動, r: 後輪駆動, 4: 4WD

# 4.
ggplot(data=mpg) + geom_point(mapping=aes(x=hwy, y=cyl))

# 5.
ggplot(data=mpg) + geom_point(mapping=aes(x=class, y=drv))
# classは車の種類、drvは駆動車輪の情報で、どちらも大小関係のある値でないから

1.2.3 グラフ作成テンプレート

# 1.
# 点が青になっていないのは、color=”blue”がaesの引数になっているから。
# color は geomの関数にする必要がある。

# 2.
# カテゴリ変数 <chr>
# 連続変数      <dbl> <int>
# カテゴリ変数 → 言葉で示される
# 連続変数     →数値で示される

# 3.
# color → だんだん濃くなる
# size  → だんだん大きくなる
# shape → Error: A continuous variable can not be mapped to shape

# 4.
ggplot(data=mpg)+geom_point(mapping=aes(x=displ, y=hwy,color=model, shape=model))
# この場合、shapeしか反映されない。カラーは赤一色
# colorとshapeの順序を入れ替えても同じ

# 5.
ggplot(data=mpg)+geom_point(mapping=aes(x=displ,
y=hwy,color=model,stroke=2.5))
# プロットサイズの大小を変更する

# 6.
# TrueとFalseで色が変わる

1.5 ファセット

# 1.
# Error in combine_vars(data, params$plot_env, vars, drop = params$drop) :
# At least one layer must contain all variables used for facetting

# 2.
# 2つの変数の組み合わせの事象が無いことを示す
ggplot(data=mpg) +
+     geom_point(mapping = aes(x = drv, y = cyl))
# とは、cyl=7のところがfacetには無い

# 3.
ggplot(data = mpg) +
geom_point(mapping = aes(x = displ, y= hwy))+
facet_grid(drv ~ .)
# ~ の前後は行と列を示す。.は何もないことを示す

1.6 幾何オブジェクト

# 1.
# 折れ線グラフ:geom_path
# 箱ひげ図:geom_boxplot()
# ヒストグラム:geom_col()
# 面グラフ:geom_area(stat=”bin”)

# 2.
ggplot(data = mpg, mapping = aes(x = displ, y = hwy, color = drv)) +
geom_point() + geom_smooth(se = FALSE)

# 3.
# 凡例をなくす
# 縦横比を一緒にしたかったから

# 5.
# 同じである

# 6.
# 左上:
ggplot(data = mpg, mapping = aes(x = displ, y = hwy)) +
geom_point() + geom_smooth(se = FALSE, show.legend = FALSE)

# 右上:
ggplot(data = mpg, mapping = aes(x = displ, y = hwy, group = drv)) +
geom_point() + geom_smooth(se = FALSE, show.legend = FALSE)

# 左中:
ggplot(data = mpg, mapping = aes(x = displ, y = hwy, color = drv)) +
geom_point() + geom_smooth(se = FALSE, show.legend = FALSE)

# 右中:
ggplot(data = mpg, mapping = aes(x = displ, y = hwy)) +
geom_point(mapping = aes(color = drv)) + geom_smooth(se = FALSE, show.legend = FALSE)

# 左下:
ggplot(data = mpg, mapping = aes(x = displ, y = hwy)) +
geom_point(mapping = aes(color = drv))

# 右下: わからない

1.7 統計変換

# 1.
# ggplot2のチートシートから、デフォルトgeomは、「geom_pointrange() 」
# 始めにデータセットを作る
data = diamonds %>%
group_by(cut) %>%
summarise(min = min(depth), max = max(depth),
median = median(depth))
ggplot(data) +
geom_pointrange(
mapping = aes(x = cut, y = median, ymin = min, ymax = max))

以下は、yminとymaxがうまくいかない
ggplot(data = diamonds) +
geom_pointrange(
mapping = aes(x = cut, y = median(depth), ymin = min(depth), ymax = max(depth)))

# 2.
# geom_colは、離散xと連続yの棒グラフ
#geom_barは、離散xと離散yの棒グラフ

# 3.
# よくわからない

# 4.
# よくわからない

# 5.
ggplot(data = diamonds) +
geom_bar(mapping = aes(x = cut, y = ..prop..))
ggplot(data = diamonds) +
geom_bar(mapping = aes(x = cut, y = ..prop.., group = 1))
ggplot(data = diamonds) +
geom_bar(
mapping = aes(x = cut, fill = color, y = ..prop..)
)

# y = ..prop.. がいけないのかな?
ggplot(data = diamonds) +
geom_bar(
mapping = aes(x = cut, fill = color)
)

1.8 位置調整

# 1.
ggplot(data = mpg, mapping = aes(x = cty, y = hwy)) +
geom_point()

# 重なってしまっていること。ジッタを入れればよい。
ggplot(data = mpg, mapping = aes(x = cty, y = hwy)) +
geom_jitter()

# 2.
# widthとheight
ggplot(data = mpg, mapping = aes(x = cty, y = hwy)) +
+     geom_jitter(width=0.4, height=0.4)

# 3.
ggplot(data = mpg, mapping = aes(x = cty, y = hwy)) +
geom_count()
# countはその点の度数に応じてプロットが大きくなる

# 4.
# わかりません

1.9 座標系

# 1.
ggplot(data = diamonds) +
geom_bar(mapping = aes(x = cut, fill = clarity))

ggplot(data = diamonds) +
geom_bar(mapping = aes(x = cut, fill = clarity)) +
coord_polar()

# 2.
ggplot(data = diamonds) +
geom_bar(mapping = aes(x = cut, fill = clarity)) +
labs( x = “New x axis label”, y = “New y axis label”,
title =”Add a title above the plot”,
subtitle = “Add a subtitle below title”)

# 3.
# coord_mapは、下行をインストールする必要がある
# install.packages(“mapproj”)
# latの表示範囲が少し違う

# 4.
ggplot(data = mpg, mapping = aes(x = cty, y = hwy)) +
geom_point() +
geom_abline() +
coord_fixed()
# 何を伝えるか:一次近似できそう
# coord_fixed() が無いと方眼にならない
# geom_abline() はx=yを表している

2章 ワークフロー:基本

2.3 関数呼び出し

# 1.
# iが違う字になっている

# 2.
library(tidyverse)

ggplot(dota = mpg) +
geom_point(mapping = aes(x=displ, y = hwy))
# dota -> data
ggplot(data = mpg) +
geom_point(mapping = aes(x=displ, y = hwy))

fliter(mpg, cyl = 8)
# fliter -> filter, = -> ==
filter(mpg, cyl == 8)

filter(diamond, carat > 3)
# diamond -> diamonds
filter(diamonds, carat > 3)

# 3.
# Help -> Keybord shortcuts help

3章 dplyrによるデータ変換

3.2.3 欠損値

# 1.
# a.
filter(flights, arr_time, sched_arr_time, arr_delay)
# これから、dep_delayが差を取っていることが分かったので、
filter(flights, arr_delay >= 60 * 2)
View(filter(flights, arr_delay >= 60 * 2))

# b.
View(filter(flights, dest == “IAH” | dest == “HOU”))

# c.
# airlinesで、対応を見ると、
# United  : UA
# American: AA
# Delta   : DL
View(filter(flights, carrier == “UA” | carrier == “AA” | carrier == “DL”))

# d.
View(filter(flights, month %in% c(7, 8, 9)))

# e.
View(filter(flights, dep_delay <= 0 & arr_delay > 120))

# f.
View(filter(flights, dep_delay > 60 & dep_delay – arr_delay > 30))

# g.
View(filter(flights, dep_time >= 000 & dep_time <=600 & arr_time >= 000 & arr_time <= 600))

# 2.
View(filter(flights, between(dep_delay, 0, 120)))

# 3.
# NAの変数は、dep?time, dep_delay, arr_time, arr_delay
# 欠航便を表す

# 4.
# そうやって組んだから
View(filter(flights, is.na(dep_time)))

3.3 arrange()で行を配置する

# 1.
arrange(df, desc(is.na(x)), x)

# 2.
View(arrange(flights, desc(arr_delay)))
View(arrange(flights, dep_time))
View(arrange(flights, desc(dep_time)))

# 3.
View(arrange(flights, distance/air_time))

# 4.
View(arrange(flights, distance))
View(arrange(flights, desc(distance)))

3.4 select()で列を選ぶ

# 1.
View(select(flights, dep_time, dep_delay, arr_time, arr_delay))

# 2.
#繰り返すと初回に出てきた時のみ表示される
View(select(flights, dep_time, arr_time, arr_delay, dep_time, dep_delay, dep_delay, arr_time, arr_delay))

# 3.
# 文字ベクトルを列挙したことと同じことになるため
# one_of(): variables in character vector
vars <- c(“year”, “month”, “day”, “dep_delay”, “arr_delay”
)
View(select(flights, one_of(vars)))

# 4.
# フィールド名にTIMEが含まれているデータを列挙する
# contains中の大文字小文字は無視するようである
View(select(flights, contains(“TIME”)))

# 大文字小文字を識別するときには、ignore.caseを用いる
View(select(flights, contains(“TIME”,ignore.case=FALSE)))
View(select(flights, contains(“time”,ignore.case=FALSE)))

3.5.1 有用な作成関数

# 1.
transmute(
flights, dep_time %/% 100 * 60 + dep_time %% 100 ,
sched_dep_time %/% 100 * 60 + sched_dep_time %% 100
)

# 2.
transmute(flights, air_time, arr_time – dep_time)

# air_time: 飛行時間[分]
# arr_time – dept_time: 時と分が混じっているので、意味をなさない

View(transmute(flights,
air_time, arr_time, dep_time,
(arr_time %/% 100 – dep_time %/% 100)*60 + (arr_time %%100 – dep_time %%100))
)
# なぜか一致しない

# 3.
View(transmute(flights,
dep_time, sched_dep_time, dep_delay))
# sched_dep_time + dep_delay = dep_time ただし、時と分が違う
View(transmute(flights,
dep_time, tmp=(sched_dep_time %/% 100)*60 + sched_dep_time %% 100 + dep_delay,
tmp %/% 60 * 100 + tmp %% 60))

# 4.
View(filter(flights, min_rank(-dep_delay)<=10))
View(arrange(filter(flights, min_rank(-dep_delay)<=10),desc(dep_delay)))

# 5.
1:3 + 1:10
# [1]  2  4  6  5  7  9  8 10 12 11
# Warning message:
# In 1:3 + 1:10 :
# longer object length is not a multiple of shorter object length

# 1:3は繰り返して適用される

# 6.
cos(x)
sin(x)
tan(x)

acos(x)
asin(x)
atan(x)
atan2(y, x)

cospi(x)
sinpi(x)
tanpi(x)

3.6.6 グループ解除

# 1.
# よくわからない

# 2.
# Q:
not_cancelled %>% count(dest)
# A:
View(group_by(not_cancelled, dest) %>%
summarize(not_cancelled = n()))

# Q:
# ‘ 以下1行にTypoあり(cancel ledと、lとlの間にスペースがある)
not_cancelled %>% count(tailnum, wt = distance)
# A:
View(group_by(not_cancelled, tailnum) %>% g
summarize(not_cancelled = sum(distance)))

# 3.
View(flights %>% filter(is.na(dep_delay)|is.na(arr_delay)))
# だと、出発したのにarr_delay=NAが含まれてしまう。(全9430件)
# 出発したのだからキャンセルではないのではないか

# 出発していないことを検出する下行ではないか
View(flights %>% filter(is.na(dep_time)))
#(全8255件)

# 4.
cancelled <- flights %>% filter(is.na(dep_time))
daily_cancelled <- group_by(cancelled, day)
count(daily_cancelled) %>% ggplot + geom_point(mapping=aes(x=day, y=n))
# なぜか8日くらいの発生が多い?

cancel_cnt <- daily_cancelled %>% summarize(cancel_cnt = n())
all_cnt <- summarize(group_by(flights, day), all_cnt = n())
dep_delay <- filter(flights, !is.na(dep_delay)) %>%
group_by(day) %>%
summarize(dep_delay = mean(dep_delay))

delay_vs_rate <- merge(
mutate(
merge(cancel_cnt, all_cnt),
cancel_rate=(cancel_cnt/all_cnt*100)
),
dep_delay
)

ggplot(data = delay_vs_rate, mapping = aes(x=dep_delay, y=cancel_rate)) + geom_point()

# 遅延が10分ほどではキャンセル率が約1.5%、20分になると約5%になる

# 5.
# ヒント
flights %>%
group_by(carrier, dest) %>%
summarize(n())

# 空港ごとに出発時の遅延と到着時の遅延の平均をとっておき、
# 個々のフライトからそれらを引く
# その後、航空会社で毎に遅延を求めることで影響分離は可能なのではないか
# やり方は分からない

mean_dep_delay <- group_by(not_cancelled, origin) %>%
summarize(m_dep_delay = mean(dep_delay))
mean_arr_delay <- group_by(not_cancelled, dest) %>%
summarize(m_arr_delay = mean(arr_delay), cnt = n())
tmp <- merge(merge(not_cancelled, mean_dep_delay),mean_arr_delay)
group_by(mutate(tmp, true_delay = arr_delay – m_dep_delay – m_dep_delay), carrier) %>%
summarize(mean(true_delay))
# F9が一番悪いということになった

# 6.
summarize(group_by(filter(not_cancelled,arr_delay<60),tailnum), cnt=n())

# 7.
# 何もしないと思うが。。。

3.7 グループごとの変更(とフィルタ)

# 1.
# 略

# 2.
delay <- group_by(not_cancelled, tailnum) %>%
summarise(delay_val = mean(dep_delay)+mean(arr_delay))
# まだ出てきてないけどorderを使う?
delay[order(delay$delay_val, decreasing = TRUE),]

# 1 N844MH        617
# 2 N911DA        562
# 3 N922EV        550

# 3.
time_vs_delay <- not_cancelled %>%
group_by(dep_time = dep_time %/% 100) %>%
summarize(arr_delay = mean(arr_delay))

ggplot(time_vs_delay) +
geom_point(mapping = aes(x = dep_time, y = arr_delay))

# 夜遅くなる方が遅延が大きくなる
# 10:00頃までに飛行すれば遅延はほとんどない

# 4.
delay_by_min <- not_cancelled %>%
mutate(arr_delay_by_min = (arr_delay %/% 100) * 60 +  arr_delay %% 100)

# 目的地ごとに計算した総遅延時間
dest_delay <- group_by(delay_by_min, dest) %>%
summarize(delay_by_dest = sum(arr_delay_by_min))

delay_data <- group_by(delay_by_min, dest, tailnum) %>%
summarize(delay_by_dest = sum(arr_delay_by_min))

dest_data <- merge(dest_delay, delay_data, by.x=”dest”, by.y=”dest”)
mutate(dest_data, rate = delay_by_dest.y / delay_by_dest.x * 100)

# 5.

lag_test <-  filter(flights, year, month, day, origin==”JFK”) %>%
mutate(lag=lag(dep_delay))
ggplot(data = lag_test) + geom_point(mapping = aes(x = dep_delay, y = lag))

# dep_delayをx軸に、lagをy軸に取ると、ほとんどは約450分を半径とした1/4円内に入る。
# 遅延が500分を超えるような便は、空港よりも飛行機に問題があると思われる
# そのため、後の便を先に出しているようである。

# 6.
# 平均と標準偏差を求める
origin_dest <- group_by(not_cancelled, origin, dest) %>%
select (origin, dest, air_time)
mean_sd <- summarise(origin_dest, mean = mean(air_time), sd=sd(air_time))
merge(origin_dest, mean_sd)
# 偏差値を求める
dev_value <- merge(origin_dest, mean_sd) %>%
mutate(dev_value = (air_time-mean)/sd * 10 + 50)

# 7.
carrier_list <- group_by(flights, dest, carrier) %>%
count() %>%
group_by(carrier) %>%
count()
sort(carrier$carrie_list)

# 9E, AA, AS, B6の順?

4章 ワークフロー:スクリプト

4.2 RStudioの診断

#1. 略
#2. 略

5章 探索的データ分析

# 5.3.3 異常値

# 1.
# xとyの分布は似ているが、zは異なる

ggplot(diamonds) +
geom_histogram(mapping = aes(x = x), binwidth = 0.5) +
coord_cartesian(xlim = c(0,10))

ggplot(diamonds) +
geom_histogram(mapping = aes(x = y), binwidth = 0.5) +
coord_cartesian(xlim = c(0,10))

ggplot(diamonds) +
geom_histogram(mapping = aes(x = z), binwidth = 0.5) +
coord_cartesian(xlim = c(0,10))

# 長さと幅は区別がつかない。高さは長さの約60%なので区別がつく

# 2.

ggplot(diamonds) +
geom_histogram(mapping = aes(x = price), binwidth = 100) +
coord_cartesian(xlim = c(0,2500))

# $1,500あたりで分布が0近くになる

# 3.
ggplot(diamonds) +
geom_histogram(mapping = aes(x = carat),binwidth=0.01)+
coord_cartesian(xlim = c(0.98,1.01))

filter(diamonds, carat==0.99) %>% count()
filter(diamonds, carat==1) %>% count()
# 0.99カラットは23個、1.00カラットは1558個
# 桁が変わるので、高級感も変わる

# 4.
ggplot(diamonds) +
geom_histogram(mapping = aes(x = z),
binwidth = 0.5)+
xlim(c(0,10))

# 特に変わりはなさそう

ggplot(diamonds) +
geom_histogram(mapping = aes(x = z), binwidth = 1)+
xlim(c(0,4))

# 棒の半分だけが表示されるということは無く、
# 一つ前のbinで切られてしまう

5.4 欠損値

# 1.
# ヒストグラム
diamonds %>%
ggplot() +
geom_histogram(mapping = aes(y)) +
xlim(c(0,10))

diamonds %>%
mutate(y = ifelse((y > 5 & y < 7) & (x < 5 | x > 7), NA, y)) %>%
ggplot() +
geom_histogram(mapping = aes(y)) +
xlim(c(0,10))
# → binにNA値以外の値が入れば、それは表示される

# 棒グラフ
diamonds %>%
mutate(cut = ifelse((cut == ‘Very Good’) & (clarity != ‘VS1’), NA, cut)) %>%
ggplot() +
geom_bar(mapping = aes(x = cut))
# → ヒストグラムと大差ない

# 2.
summarize(diamonds, mean(x))
summarize(diamonds, mean(x, na.rm = TRUE))
summarize(diamonds, mean(x, na.rm = FALSE))
summarize(diamonds, sum(x))
summarize(diamonds, sum(x, na.rm = TRUE))
summarize(diamonds, sum(x, na.rm = FALSE))
# → どちらも特に変わらないようである

5.5.1 カテゴリ変数と連続変数

# 1.
mutate(flights, flg = is.na(dep_time)) %>%
ggplot(mapping = aes(x = sched_dep_time)) +
geom_freqpoly(mapping = aes(color = flg),binwidth = 100)

# 2.
ggplot(diamonds) +
geom_point(mapping = aes(x = carat, y = price))
ggplot(diamonds) +
geom_boxplot(mapping = aes(x = cut, y = price))
ggplot(diamonds) +
geom_boxplot(mapping = aes(x = color, y = price))
ggplot(diamonds) +
geom_boxplot(mapping = aes(x = clarity, y = price))
ggplot(diamonds) +
geom_point(mapping = aes(x = depth, y = price))
ggplot(diamonds) +
geom_point(mapping = aes(x = table, y = price))
ggplot(diamonds) +
geom_point(mapping = aes(x = x, y = price))
ggplot(diamonds) +
geom_point(mapping = aes(x = y, y = price))
ggplot(diamonds) +
geom_point(mapping = aes(x = z, y = price))
# 価格を予想するのにどの変数が重要かはよく分からない
# 仮にcaratが重要だとする
ggplot(diamonds) +
geom_boxplot(mapping = aes(x = cut, y = carat))
# 一番品質の低いFairの50百分位が最も高い位置にある

# 3.
install.packages(“ggstance”)
library(ggstance)
ggplot(data = mpg) +
geom_boxploth(
mapping = aes(
x = hwy,
y = reorder(class, hwy, FUN = median)
)
)
# 違いが全く分からない

# 4.
install.packages(“lvplot”)
library(lvplot)
ggplot(data = diamonds, mapping = aes(x = cut, y = price)) +
geom_lv()
# 横幅で度数の比率を表している
# Fairには高額の度数比率が少ない
# Idealは、高額の度数比率が少なくならない

# 5.
ggplot(data = diamonds, mapping = aes(x = cut, y = price)) +
geom_violin()
# geom_violin について
# 利点:カテゴリ毎の分布の割合が分かる
# 欠点:実数量は分からない

ggplot(data = diamonds, mapping = aes(x = price, fill=cut)) +
geom_histogram()
# 積層geom_histgramについて
# 利点:bin毎のトータルの実数量が分かる。bin内のカテゴリの割合も分かる
# 欠点:一つのカテゴリ毎分布は分かりにくい

ggplot(data = diamonds, mapping = aes(x = price)) +
geom_freqpoly(mapping = aes(color = cut), binwidth = 500)
# 色つきのgeom_freqpolyについて
# 利点:カテゴリ毎にbin毎の数量が分かる
# 欠点:線が混み入ってくると分かりにくくなる

# 5.
#install.packages(“ggbeeswarm”)
library(ggbeeswarm)
View(available.packages())
help(package=”ggbeeswarm”)

5.5.2 2つのカテゴリ変数

# 1.
diamonds %>%
count(color, cut) %>%
ggplot(mapping = aes(x = color, y = cut)) +
geom_tile(mapping = aes(fill = n > 2000))

# 2.
not_cancelled <- flights %>%
filter(!is.na(dep_delay), !is.na(arr_delay))

flights %>%
group_by(dest, month) %>%
summarize(mean = mean(arr_delay)) %>%
ggplot(mapping = aes(x = month , y =dest)) +
geom_tile(mapping = aes(fill = mean))
# 目的地が多すぎてわからない
# 目的地を絞る、何かしらソートするなど

# 3.
diamonds %>%
count(color, cut) %>%
ggplot(mapping = aes(x = cut, y = color)) +
geom_tile(mapping = aes(fill = n))
# 視線は左右に移しやすく、後者の変化方向が左右だから

5.5.3 2つの連続変数

# 1.
# 何がなんだかわからない
ggplot(data = smaller, mapping = aes(x = carat, y = price)) +
geom_area(mapping = aes(group = cut_width(carat, 0.1)))

ggplot(data = smaller, mapping = aes(x = price, y = carat)) +
geom_area(mapping = aes(group = cut_width(price, 1000)))

ggplot(data = smaller, mapping = aes(x = carat, y = price)) +
geom_area(mapping = aes(group = cut_number(carat, 20)))

ggplot(data = smaller, mapping = aes(x = price, y = carat)) +
geom_area(mapping = aes(group = cut_number(price, 10)))

# 2.
ggplot(data = smaller, mapping = aes(x = price, y = carat)) +
geom_boxplot(mapping = aes(group = cut_number(price, 100)))

# 3.
ggplot(data = smaller, mapping = aes(x = carat, y = price)) +
geom_boxplot(mapping = aes(group = cut_number(carat, 20)))
# 小さなダイヤには高額のものが点在する
# 大きなダイヤには定額のものが点在する

# 4.
# 問の意味がよく分からない
ggplot(data = diamonds) +
geom_point(mapping = aes(x = cut, y = price),
alpha = 1/100)

# 5.
ggplot(data = diamonds) +
geom_point(mapping = aes(x = x, y = y)) +
coord_cartesian(xlim = c(4,11), ylim = c(4,11))

ggplot(data = diamonds, mapping = aes(x = x, y = y)) +
geom_boxplot(mapping = aes(group = cut_number(x, 10)))+
coord_cartesian(xlim = c(4,11), ylim = c(4,11))
# binの境目の情報がどちらかに丸められてしまうためであろう
# bin数を多くすれば散布図に近づく

7章 tibbleのtibble

7.4 古いコードとの関わり

# 1.
print(mtcars)
print(tb)
# printした時に A tibble: と示される

# 2.
df <- data.frame(abc = 1, xyz = “a”)
df$x
df[,”xyz”]
df[,c(“abc”,”xyz”)]

tb <- tibble(abc = 1, xyz = “a”)
tb$x
tb$xyz
tb[,”xyz”]
tb[,c(“abc”,”xyz”)]
# data.frameは、コラム名を勝手に補完してしまう

# 3.
var <- “mpg”
# 問の意味が分からない

# 4.
# データフレームと書いてあるが、tibbleのことだと思う
# a.
annoying <- tibble(`1` = 1:10, `2` = `1` * 2 + rnorm(length(`1`)))
# b.
ggplot() +
geom_point(data = annoying, mapping = aes(x = `1`, y = `2`))
# c.
tb = mutate(annoying, `3` = `2` / `1`)
# d.
names(tb)
names(tb)[1] <- “one”
names(tb)[2] <- “two”
names(tb)[3] <- “three”
names(tb)

# 5.
?tibble::enframe()
tb <- tribble(
~x, ~y, ~z,
#–|–|—-
“a”, 2, 3.6,
“b”, 1, 8.5,
“c”, 5, 2.8
)
df <- tibble::enframe(tb)
# tibble::enframeはTibbleをタイトルを1列目、
# データベクトルを2列目とするテーブルに変換する

# 6.
?print.tbl_df
# この結果より、n_extraであることがわかる
library(nycflights13)
print(flights, n_extra = 3)

8章 readrによるデータインポート

8.2.1 基本Rとの比較

# 1.
read_delim(“../r4ds-master/data/|.csv”,”|”)

# 2.
?read_csv
# col_names, col_types, locale, na, quoted_na, quote, trim_ws, n_max, guess_max, protgress

# 3.
?read_fwf
# col_positionsか

# 4.
read_delim(“x,y\n,’a,b'”,”,”,quote=”\'”)

# 5.
read_csv(“a,b\n1,2,3\n4,5,6”)
# → タイトルが2つなのにデータが3つずつある
# → 3つ目のデータは無視される

read_csv(“a,b,c\n1,2\n1,2,3,4”)
# → タイトルが3つなのに、データが2つまたは4つである
# → 2つのデータは3つ目にNAが入る4つ目のデータは無視される

read_csv(“a,b\n\”1”)
# → 1が数値として認識される。2データ目はNA

read_csv(“a,b\n1,2\na,b”)
# → chrとして認識されるだけである

8.3.4 日付、日付時刻、時刻

# 1.
?locale
# date_names

# 2.
parse_number(“$123,456,789.0”,
locale = locale(decimal_mark = “,”, grouping_mark = “,”))
# Error: `decimal_mark` and `grouping_mark` must be different

parse_number(“$123,456.7890”,
locale = locale(decimal_mark = “,”))
# grouping_markのデフォルト値は不明だが、カンマは無視される

parse_number(“$123,456,789.0”,
locale = locale(grouping_mark = “.”))
# カンマが小数点になる

# 3.
# date_format: デフォルトの日付の形式
# time_format: デフォルトの時刻の形式

parse_date(“01/02/2013”, locale = locale(date_format = “%d/%m/%Y”))
parse_time(“08:40:23”, locale = locale(time_format = “%S:%M:%H”))

# 4.
# 分からない

# 5.
# read_csv2はセパレータとしてセミコロンを用いる

# 6.
# ISO-8859-1か?

# 7.
d1 <- “January 1, 2010”
d2 <- “2015-Mar-07”
d3 <- “06-Jun-2017”
d4 <- c(“August 15 (2015)”, “July 1 (2015)”)
d5 <- “12/30/14” # Dec 30, 2014
t1 <- “1705”
t2 <- “11:15:10.12 PM”

parse_date(d1,”%B %d, %Y”)
parse_date(d2,”%Y-%b-%d”)
parse_date(d3,”%d-%b-%Y”)
parse_date(d4,”%B %d (%Y)”)
parse_date(d5,”%m/%d/%y”)
parse_time(t1,”%H%M”)
parse_time(t2,”%I:%M:%OS %p”)

9章 tidyrによるデータ整理

9.2 整理データ

# 1.
# table1は、1つの列に、国名、年、症例数、人口が記されている
# Table2は、type列により、症例行か人口行かに区別されている
# Table3は、人口当たりの症例発生率をスラッシュを使い、1列で表現している
# Table4は、aに症例数、bに人口を記載している

# 2.
# table2
table2a <- filter(table2, type == “cases”)
table2b <- filter(table2, type == “population”)
inner_join(table2a, table2b, by = c(“country”, “year”)) %>%
mutate(rate = count.x/count.y * 10000)

# table3
table3a <- as_tibble(t(sapply(strsplit(table3$rate,”/”),
function(x) x[])))
table3b <- cbind(table3, table3a)
table3c <-
type_convert(table3b,
cols(V1 = col_double(), V2 = col_double()))
mutate(table3c, rate1 = V1 / V2 * 10000)

# table4
table4c <- inner_join(select(table4a,c(country,cases = `1999`)),
select(table4b,c(country,population = `1999`)),
by = “country”) %>%
mutate(year = 1999)
table4d <- inner_join(select(table4a,c(country,cases = `2000`)),
select(table4b,c(country,population = `2000`)),
by = “country”) %>%
mutate(year = 2000)
table4e <- rbind(table4c, table4d)
table4e[order(table4e$country),] %>%
mutate(rate = cases/population * 10000)

# 所要行数の多い問題を難しいとすると、最も難しいのはtable4
# 名前の付け替えや行列の分離・結合があるため

# 3.
ggplot(type_convert(filter(table2, type == “cases”)), aes(year, count)) +
geom_line(aes(group = country), color = “gray50”) +
geom_point(aes(color = country))

9.3 広げたり集めたり

# 1.
stocks <- tibble(
year   = c(2015, 2015, 2016, 2016),
half   = c(   1,    2,    1,    2),
return = c(1.88, 0.59, 0.92, 0.17)
)
stocks
stocks %>%
spread(year, return) %>%
gather(“year”, “return”, `2015`:’2016′)
# 中間段階で、西暦がタイトルになるため、dblからChrに変換されてしまう

# convert = TRUEオプションにより、spread後に適切な型になる
stocks <- tibble(
year   = c(2015, 2015, 2016, 2016),
half   = c(   1,    2,    1,    2),
return = c(1.88, 0.59, 0.92, 0.17)
)
stocks
stocks %>%
spread(year, return) %>%
gather(“year”, “return”, `2015`:’2016′, convert=TRUE)

# yearがdblからintにはなったが、少なくともchrではなくなった

# 2.
table4a %>%
gather(1999, 2000, key = “year”, value = “cases”)
# 1999と2000が非構文的名前なので、バッククォートで括る必要がある

# 3.
people <- tribble(
~name,            ~key,     ~value,
#—————–|———|——
“Phillip Woods”,   “age”,     45,
“Phillip Woods”,   “height”, 186,
“Phillip Woods”,   “age”,     50,
“Jessica Cordero”, “age”,     37,
“Jessica Cordero”, “height”, 156
)
people %>% spread(key, value)
# Error: Duplicate identifiers for rows (1, 3)

people2 <- tribble(
~name,            ~key,     ~value, ~cnt,
#—————–|———|——|——
“Phillip Woods”,   “age”,     45,   1,
“Phillip Woods”,   “height”, 186,   1,
“Phillip Woods”,   “age”,     50,   2,
“Jessica Cordero”, “age”,     37,   1,
“Jessica Cordero”, “height”, 156,   1
)
people2 %>% spread(key, value)

# 4.
preg <- tribble(
~pregnant, ~male, ~female,
“yes”,    NA,      10,
“no”,    20,      12
)
preg %>% gather(male, female, key = “M/F”, value = “count”)

9.4.2 結合

# 1.
tibble(x = c(“a,b,c”, “d,e,f,g”, “h,i,j”)) %>%
separate(x, c(“one”, “two”, “three”))
tibble(x = c(“a,b,c”, “d,e”, “f,g,i”)) %>%
separate(x, c(“one”, “two”, “three”))
# —
tibble(x = c(“a,b,c”, “d,e,f,g”, “h,i,j”)) %>%
separate(x, c(“one”, “two”, “three”), extra = “merge”)
# extra = “merge” の場合、引数が無視されないで無理に使われる
tibble(x = c(“a,b,c”, “d,e”, “f,g,i”)) %>%
separate(x, c(“one”, “two”, “three”), fill = “left”)
# fill = “left”の場合、NAがあったときに左詰めになる

# 2.
table5 %>%
unite(new, century, year, sep = “”, remove = FALSE)
table5 %>%
unite(new, century, year, sep = “”, remove = TRUE)
table3 %>%
separate(rate, into = c(“cases”, “population”), remove = FALSE)
table3 %>%
separate(rate, into = c(“cases”, “population”), remove = TRUE)
# remove = TRUEにすると、元の列が消される

# 3.
# extractは正規化表現でグルーピングして分割する
# 結合云々については分からない
table3 %>%
separate(rate, into = c(“cases”, “population”))
table3 %>%
extract(rate, into = c(“cases”, “population”), regex = ‘(.+)/(.+)’)

9.5 欠損値

# 1.
# fill(data, …, .direction = c(“down”, “up”))
# spread(data, key, value, fill = NA, convert = FALSE, drop = TRUE, sep = NULL)
# complete(data, …, fill = list())

# 2.
treatment %>%
fill(person, .direction = “up”)
# “up”の場合、下に書かれている値で補完される

9.6 ケーススタディ

# 1.
# 症例数を数えるということではna.rm = TRUEは妥当である
# 値のあるべきセルにNAがあるので、明示的欠損値である(9.5節より)
# 0は調査した結果、その症例が無かったことを示し、NAは調査しなかったかもしれない

# 2.
# Warning message:
# 次のseparate文でエラーが発生する
# Expected 3 pieces. Missing pieces filled with `NA` in 2580 rows [73467, 73468, 73469, 73470, 73471, 73472, 73473, 73474, 73475, 73476, 73477, 73478, 73479, 73480, 73481, 73482, 73483, 73484, 73485, 73486, …].

# 3.
select(who, country, iso2, iso3) %>%
mutate(paste = paste(paste(country, iso2), iso3)) %>%
select(paste) %>%
count(paste)
# これで、重複が無いことを確認する

# 4.
# 各国
ggplot(who5) +
geom_bar(mapping = aes(x=country, y=cases), stat=”identity”)
ggplot(who5) +
geom_bar(mapping = aes(x=age, y=cases), stat=”identity”)
ggplot(who5) +
geom_bar(mapping = aes(x=sex, y=cases), stat=”identity”)

10章 dplyrによる関係データ

10.2 nycflights13

# 1.
# ・airportのname, lat, lon
# ・planesのtype
# ・weatherのorigin, year, month, day, hour, temp, dewp, humid

# 2.
# weatherのoriginが半白丸、airportsのfaaが矢印で結びつく

# 3.
# arr_year, arr_month, arr_day

# 4.
# data frame: special_day
# composite primary key: date + faa
# field: count

10.3 キー

# 1.
lfights_ak <- mutate(flights, row_number())

# 2.
# a.
Lahman::Batting %>%
count(playerID, yearID, stint) %>%
filter(n > 1)

# b.
# install.packages(“babynames”)
# library(babynames)
babynames::babynames %>%
count(year, sex, name) %>%
filter(nn > 1)

# c.
# install.packages(“nasaweather”)
# library(nasaweather)
nasaweather::atmos %>%
count(lat,long, year, month) %>%
filter(n > 1)

# d.
# install.packages(“fueleconomy”)
# library(fueleconomy)
fueleconomy::vehicles %>%
count(id) %>%
filter(n > 1)

# e
ggplot2::diamonds %>%
count(carat, cut, color, clarity, depth, table, price, x, y, z) %>%
filter(n > 1)
# キーはない

# 3.
# 略

10.4.5 キーの列を定義する

# 1.
airports %>%
semi_join(flights, c(“faa” = “dest”)) %>%
ggplot(aes(lon, lat)) +
borders(“state”) +
geom_point() +
coord_quickmap()

not_cancelled <- flights %>%
filter(!is.na(dep_delay), !is.na(arr_delay))
mean_arr_delay <- group_by(not_cancelled, dest) %>%
summarize(m_arr_delay = mean(arr_delay))

airports %>%
inner_join(mean_arr_delay , c(“faa” = “dest”)) %>%
ggplot(aes(lon, lat)) +
borders(“state”) +
geom_point(aes(color = m_arr_delay)) +
coord_quickmap()

# 2.
flights_loc <- flights %>%
left_join(airports, c(“origin” = “faa”)) %>%
left_join(airports, c(“dest” = “faa”))
View(flights_loc)

# 3.
mean_arr_delay
planes
mean_delay <- group_by(not_cancelled, tailnum) %>%
summarise(mean_delay = mean(arr_delay))
inner_join(mean_delay, planes) %>%
mutate(age = 2018-year) %>%
select(age, mean_delay) %>%
ggplot() +
geom_point(aes(age, mean_delay))
# 10~30年くらいの機体に遅延が多い事例が散見される

# 4.
#flightsとwetherのjoinを取り、遅延と温度等の相関をとればよい
#組合せ数が膨大になるので省略

# 5.
not_cancelled %>%
filter(year == 2013, month == 06, day == 13) %>%
ggplot() +
geom_point(aes(dep_time, arr_delay))
# 竜巻の影響で、時刻が遅くなるにつれて遅延する飛行機の数や遅延時間が増えた
# 天候との相互参照は省略

10.5 フィルタジョイン

# 1.
View(filter(flights, is.na(tailnum)))
# tailnumが欠損しているフライトは、飛行記録が無いので、
# 機体手配ができなかったのではないか

View(flights %>%
anti_join(planes, by = “tailnum”) %>%
count(tailnum, sort = TRUE))
# planesにマッチするレコードがない機体記号に共通するのは、
# 下2桁がMQまたはAAである
View(flights %>%
semi_join(flights %>%
anti_join(planes, by = “tailnum”) %>%
count(tailnum, sort = TRUE)%>%
filter(!is.na(tailnum))))
# carrierが機体を登録していない?

# 2.
flights %>%
semi_join(flights %>%
count(tailnum) %>%
filter(!is.na(tailnum), n>100))

# 3.
fe <- inner_join(fueleconomy::vehicles, fueleconomy::common)
inner_join(head(count(fe,model,sort=TRUE),1), fe)

# 4.
not_cancelled <- flights %>%
filter(!is.na(dep_delay), !is.na(arr_delay))
arr_delay_total <- group_by(not_cancelled, year, month, day) %>%
summarise_each(funs(sum), arr_delay)
arr_delay_total[order(arr_delay_total$arr_delay,decreasing = TRUE),]
# 2013/03/08が最も遅延の和が大きい(68518)
weather %>%
filter(year==2013, month==3, day==8) %>%
ggplot() +
geom_point(mapping = aes(x=hour, y=temp))
# 昼間が温度が低く(約31?)、夜間が高い(42?23時)、
# でも、他の日とあまり変わらない
weather %>%
filter(year==2013, month==3, day==8) %>%
ggplot() +
geom_point(mapping = aes(x=hour, y=humid))
# 湿度はずっと90%くらい。高い
weather %>%
filter(year==2013, month==3, day==8) %>%
ggplot() +
geom_point(mapping = aes(x=hour, y=wind_dir))
# 風の方向はずっと北風
weather %>%
filter(year==2013, month==3, day==8) %>%
ggplot() +
geom_point(mapping = aes(x=hour, y=wind_speed))
# 風速は18mphくらい。そんなに速くない
weather %>%
filter(year==2013, month==3, day==8) %>%
ggplot() +
geom_point(mapping = aes(x=hour, y=wind_gust))
# 突風も20mphくらい。そんなに速くない
# よって、湿度が高いことが原因のようである

# 5.
anti_join(flights, airports, by = c(“dest”=”faa”)) %>%
count(dest)
# BQN, QSE, SJU, STTがairportsに登録されておらず、
# 一日に20便くらいはそれらの空港に飛んでいる

anti_join(airports, flights, by = c(“faa”=”dest”))
# 1300もの空港に対する飛行データがflightsでは欠損している

# 6.
flights %>%
count(carrier, tailnum) %>%
count(tailnum) %>%
filter(nn > 1) %>%
inner_join(count(flights,carrier, tailnum))%>%
spread(key = carrier, value = n)
# 1つの機体を複数の航空会社で運航している例がある
# 9EとEVがN146PQを含む8機体、、DLとFLがN933A%を含む9機体を共有している

11章 stringrによる文字列

11.2.4 ロケール

# 1.

paste(“1st”, “2nd”, “3rd”)
paste0(“1st”, “2nd”, “3rd”)

# paste は、デフォルトで半角スペースで文字列をつなぐ
# paste0は、文字と文字の間に何も入れないで文字列をつなぐ

# stringr の関数はstr_c
str_c(“1st”, “2nd”, “3rd”, sep=’,’)

# 2.
# sep は普通の区切り文字
# collapse はベクトルの要素に対して使う?
paste0(c(“1st”, “2nd”, “3rd”), collapse = “, “)

# 3.
# 偶数のときに真ん中2文字を取る場合
strData <- “abcde”
str_sub(strData, (str_length(strData)+1)%/%2, (str_length(strData)+2)%/%2)
strData <- “abcdef”
str_sub(strData, (str_length(strData)+1)%/%2, (str_length(strData)+2)%/%2)

#偶数のときに空白にする場合
strData <- “abcde”
str_sub(strData, (str_length(strData)+1)%/%2,
(str_length(strData)+1 – ((str_length(strData)+1)%%2)*2)%/%2)
strData <- “abcdef”
str_sub(strData, (str_length(strData)+1)%/%2,
(str_length(strData)+1 – ((str_length(strData)+1)%%2)*2)%/%2)

# 4.
# str_wrapは、widthで定めた文字数ごとに改行コードを入れ、
# さらに、先頭にindentで指定した文字数、2行目以降の先頭にexdentで
# 指定した文字数の空白文字列を入れる。ただし、単語の途中で改行はされない
str_wrap(string = “abcdefghij klmn aqk;lkiawi ai”, width = 8,
indent = 2, exdent = 3)

# 5.
# str_trim は文字列前後の空白を除去する
str_trim(”   abc “)
# str_trimの逆操作はstr_pad
str_pad(“abc”,30,”both”)

# 6.
c(“a”) %>%
str_c(collapse = “, “)
c(“a”, “b”) %>%
str_c(collapse = “, “)
c(“a”, “b”, “c”) %>%
str_c(collapse = “, “)

11.3.1 基本マッチ

# 1.
# “\”  だと、エスケープ文字であって、その対象がない
# “\\” だと、正規表現で\であり、エスケープを示すが、その対象がない
# “\\\”だと、正規表現の\とエスケープ文字を意味するが、エスケープ対象がない

# 2.
str_view(“\”‘\\”,”\”‘\\\\”)

# 3.
# “.?.?.?” (ただし?は任意の1文字)のような文字とマッチする
str_view(“.a.1.B”,”\\..\\..\\..”)

11.3.2 アンカー

# 1.
str_view(“$^$”,”\\$\\^\\$”)

# 2.
# a.
stringr::words[grep(“^y”, stringr::words)]
# b.
stringr::words[grep(“x$”, stringr::words)]
# c.
stringr::words[grep(“^…$”, stringr::words)]
# d.
stringr::words[grep(“…….”, stringr::words)]

11.3.3 文字のクラスと候補

# 1.
# a.
stringr::words[grep(“^(a|i|u|e|o)”,stringr::words)]
# b.
stringr::words[grep(“^[^(a|i|u|e|o)]”,stringr::words)]
# c.
stringr::words[grep(“[^e]ed$”,stringr::words)]
# d.
stringr::words[grep(“(ing|ize)$”,stringr::words)]

# 2.
stringr::words[grep(“ie”,stringr::words)]
stringr::words[grep(“cie”,stringr::words)]
# “science”と”society”の例外がある

# 3.
# 続く
stringr::words[grep(“q[^u]”,stringr::words)]
stringr::words[grep(“qu”,stringr::words)]

# 4.
stringr::words[grep(“re$”,stringr::words)]
stringr::words[grep(“our$”,stringr::words)]
stringr::words[grep(“ise$”,stringr::words)]
stringr::words[grep(“yse$”,stringr::words)]
stringr::words[grep(“ence$”,stringr::words)]
stringr::words[grep(“ogue$”,stringr::words)]

# 5.
str_view(“080-1234-5678″,”\\d\\d\\d-\\d\\d\\d\\d-\\d\\d\\d\\d”)

11.3.4 繰り返し

# 1.
str_view(x,”C?”)
str_view(x,”C{0,1}”)

str_view(x,”C+”)
str_view(x,”C{1,}”)

str_view(x,”C*”)
str_view(x,”C”)

# 2.
# a.: 1文字以上の任意の文字列
str_view(“abc”,”^.*$”)

# b.: 中括弧で囲まれた1文字以上の任意の文字列
str_view(“{abc}”,”\\{.+\\}”)

# c.: 4桁、2桁、2桁の数字がハイフンでつながれた文字列
str_view(“0123-45-67″,”\\d{4}-\\d{2}-\\d{2}”)

# d.: バックスラッシュ4本からなる文字列
str_view(“\\\\\\\\”,”\\\\{4}”)

# 3.
# a.
stringr::words[grep(“^[^a^i^u^e^o][^a^i^u^e^o][^a^i^u^e^o]”,stringr::words)]

# b.
stringr::words[grep(“(a|i|u|e|o)(a|i|u|e|o)(a|i|u|e|o)”,stringr::words)]

# c.
stringr::words[grep(“(a|i|u|e|o)[^a^i^u^e^o](a|i|u|e|o)[^a^i^u^e^o]”,stringr::words)]
str_view(stringr::words,
“(a|i|u|e|o)[^a^i^u^e^o](a|i|u|e|o)[^a^i^u^e^o]”,
match = TRUE)

# 4.
# このサイトは安全に接続できませんだそうである

11.3.5 グループ化と後方参照

# 1.
# a.: 任意の1文字の3連続
str_view(“abcaabbccaaaabbbccc”, “(.)\\1\\1”)

# b.: abbaのようにn文字目とn+3文字目、n+1文字目とn+2文字目が同じ文字列
str_view(“xabbax”, “(.)(.)\\2\\1”)

# c.: ababのようにn文字目とn+2文字目、n+1文字目とn3文字目が同じ文字列
str_view(“xababx”, “(..)\\1”)

# d.: 任意の文字がn、n+2、n+4文字目に現れる文字列
str_view(“axbxcxd”, “(.).\\1.\\1”)

# e.: ある3文字の並びが任意の文字数だけ離れて逆順に現れる文字列
str_view(“xabcdefhijcbay”, “(.)(.)(.).*\\3\\2\\1”)

# 2.
# a.
str_view(“abcdefga”, “^(.).*\\1$”)

# b.
str_view(“church”, “(.*).*\\1”)

# c.
str_view(“eleven”, “(.).*\\1.*\\1”)

11.4.1 マッチの可否

# 1.
# a.
str_view(words,”^x|x$”, match = TRUE)

a <- str_detect(words, “^x”)
b <- str_detect(words, “x$”)
words[a | b]

# b.
str_view(words, “^[aeiou].*[^aeiou]$”, match = TRUE)

a <- str_detect(words, “^[aeiou]”)
b <- str_detect(words, “[^aeiou]$”)
words[a & b]

# c.
str_view(words, paste0(“^”,
“(?=.*(a.*e|e.*a))(?=.*(a.*i|i.*a))(?=.*(a.*o|o.*a))”,
“(?=.*(a.*u|u.*a))(?=.*(e.*i|i.*e))(?=.*(e.*o|o.*e))”,
“(?=.*(e.*u|u.*e))(?=.*(i.*o|o.*i))(?=.*(i.*u|u.*i))”,
“(?=.*(o.*u|u.*o))”), match = TRUE)

a     <- words[grep(“a”, words)]
ae    <- a[grep(“e”, a)]
aei   <- ae[grep(“i”, ae)]
aeio  <- aei[grep(“o”, aei)]
aeiou <- aeio[grep(“u”, aeio)]
aeiou

# d.
# 母音数が一番多い単語
sort(str_count(words, “a”) +
str_count(words, “e”) +
str_count(words, “i”) +
str_count(words, “o”) +
str_count(words, “u”),
decreasing = TRUE)

words[order(str_count(words, “a”) +
str_count(words, “e”) +
str_count(words, “i”) +
str_count(words, “o”) +
str_count(words, “u”),
decreasing = TRUE)][1]

11.4.2 マッチの抽出

# 1.
colors <- c(
“Red”, ” red”, “Orange”, ” orange”, “Yellow”, ” yellow”,
“Green”, ” green”, “Blue”, ” blue”, “Purple”, ” purple”
)
color_match <- str_c(colors, collapse = “|”)
str_view_all(str_subset(sentences, color_match), color_match)

# 2.
# a.
str_extract(sentences, “^[0-9A-Za-z]+ “)

# b.
str_extract(sentences[str_count(sentences, “.*ing”) >= 1],
“^[0-9A-Za-z]+ing| [0-9A-Za-z]+ing”)

# c.
# 複数を表すsと単語の最後のsが分からない。childとchildren等も抽出できない
str_extract(sentences[str_count(sentences, “[0-9A-Za-z]+s | [0-9A-Za-z]+s\\.”) >= 1],
“[0-9A-Za-z]+s | [0-9A-Za-z]+s\\.”)

11.4.3 グループのマッチ

# 1.
num <- “(number) ([^ ]+)”
has_num <- sentences %>%
str_subset(num) %>%
head(10)
has_num
# sentencesの中にはなさそう

# 2.
apos <- “([^ ]+)'([^ ]+)”
has_apos <- sentences %>%
str_subset(apos) %>%
head(10)
has_apos %>%
str_match(apos)

11.4.4 置換マッチ

# 1.
str_replace_all(“a/b/c”, “/”, “\\\\”)

# 2.
str_replace_all(sentences, c(
“A” = “a”,
“B” = “b”,
“C” = “c”,
“D” = “d”,
“E” = “e”,
“F” = “f”,
“G” = “g”,
“H” = “h”,
“I” = “i”,
“J” = “j”,
“K” = “k”,
“L” = “l”,
“M” = “m”,
“N” = “n”,
“O” = “o”,
“P” = “p”,
“Q” = “q”,
“R” = “r”,
“S” = “s”,
“T” = “t”,
“U” = “u”,
“V” = “v”,
“W” = “w”,
“X” = “x”,
“Y” = “y”,
“Z” = “z”
)) %>%
head(10)

# 3.
str_replace(words, “(^.)(.*)(.$)”,”\\3\\2\\1″)
replace_ok <-
words[match(words,str_replace(words, “(^.)(.*)(.$)”,”\\3\\2\\1″))]
replace_ok[!is.na(replace_ok)]

11.4.5 分割

# 1.
x <- “apples, pears, and bananas”
str_split(x, ” “)[[1]]
str_split(x, boundary(“word”))[[1]]

# 2.
# wordsで分割すると、カンマが除去される

# 3.
str_split(x, “”)[[1]]
# 1文字ずつ分割される

11.5 他の種類のパターン

# 1.
str_view(“\\\\”, “.*”)
# fixedではわからない

# 2.
df <- paste(sentences, collapse = ” “) %>%
str_extract_all(boundary(“word”)) %>%
unlist() %>%
tibble(name = .) %>%
group_by(name) %>%
count()
df[order(df$n, decreasing=TRUE),] %>%
head(5)

11.7 stringi

# 1.
# library(stringi)
# a.
apropos(“stri_count”)
stri_count_words(str_extract_all(sentences[1], boundary(“word”)))

# b.
apropos(“stri_du”)
?stri_duplicated
df <- paste(sentences, collapse = ” “) %>%
str_extract_all(boundary(“word”)) %>%
unlist()
df[stri_duplicated(df)]

# c.
apropos(“random”)
apropos(“rand”)
stri_rand_shuffle
stri_rand_shuffle(sentences[[1]])

# 2.
?stri_sort
# locale= で制御する?

12章 forcatsでファクタ

12.3 総合的社会調査

# 1.
df <- gss_cat %>%
count(rincome)
ggplot(data = gss_cat) +
geom_bar(aes(rincome)) +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
# 理解しがたいのはX軸ラベルが重なるから?金額順に並んでいないため?
income_levels <- c(
“Lt $1000″,”$1000 to 2999″,”$3000 to 3999″,”$4000 to 4999”,
“$5000 to 5999″,”$6000 to 6999″,”$7000 to 7999”, “$8000 to 9999”,
“$10000 – 14999″,”$15000 – 19999″,”$20000 – 2499″,”$25000 or more”,
“Refused”,”Don’t know”,”No answer”,”Not applicable”
)
gss_cat2 <- transform(gss_cat, rincome=factor(rincome, levels = income_levels))
ggplot(data = gss_cat2) +
geom_bar(aes(rincome)) +
theme(axis.text.x = element_text(angle = 90, hjust = 1))

# 2.
gss_cat %>% count(relig, sort = TRUE)
# Protestant
gss_cat %>% count(partyid, sort = TRUE)
# Independent

# 3.
gss_cat %>% count(relig, denom) %>%
ggplot() +
geom_point(mapping = aes(x = relig, y = denom)) +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
# 宗派はProtestantに適用される

12.4 ファクタ順序の変更

# 1.
ggplot(gss_cat) +
geom_bar(aes(x = tvhours))
# 度数が少ないので有用だと思う

# 2.
# 型がfctのファクタを調査する
gss_cat
# fctは、marital, race, rincome, partyid, relig, denomである
levels(gss_cat$marital)
# maritalは定まっているようである
levels(gss_cat$race)
# raceも定まっているようである
levels(gss_cat$rincome)
# rincomeも定まっている
levels(gss_cat$partyid)
# partyidも定まっている
levels(gss_cat$denom)
# demonも定まっている

# 3.
# Y軸は、順序が先のものが下に来るから

12.5 ファクタ水準の変更

# 1.
df <- gss_cat %>% group_by(year, denom) %>% count()
ggplot(df) +
geom_line(aes(year, n, group = denom, color=denom))

# 2.
gss_cat %>% group_by(rincome) %>% count()
gss_cat %>% mutate(rincome = fct_collapse(rincome,
other = c(“No answer”, “Don’t know”, “Refused”, “Not applicable”),
“$25000 or more” = c(“$25000 or more”),
“$1000 to $24999” = c(“$20000 – 24999″,”$15000 – 19999”,
“$10000 – 14999″,”$8000 to 9999″,”$7000 to 7999”,
“$6000 to 6999″,”$5000 to 5999″,”$4000 to 4999”,
“$3000 to 3999″,”$1000 to 2999”),
“Lt $1000” = c(“Lt $1000”)
)) %>% group_by(rincome) %>% count()

13章 lubridateによる日付と時刻

13.2.3 他の型から作成

# 1.
ymd(c(“2010-10-10”, “bananas”))
# NAになる

# 2.
?today
today(“Australia/Canberra”)
# 時差があるから

# 3.
d1 <- “January 1, 2010”
d2 <- “2015-Mar-07”
d3 <- “06-Jun-2017”
d4 <- c(“August 19 (2015)”, “July 1 (2015)”)
d5 <- “12/30/14” # Dec 30, 2014

?lubridate

strptime(d1, format = “%B %d, %Y”, tz = “America/New_York”)
mdy(d1)
parse_date(d1,”%B %d, %Y”)
# lubridateによるd1の変換の仕方がわからない

ymd(d2)

dmy(d3)

strptime(d4, “%B %d (%Y)”)
parse_date(d4, “%B %d (%Y)”)
# lubridateによるd4の変換の仕方がわからない

mdy(d5)

13.3.3 要素を設定する

# 1.
# 1月
flights_dt %>% filter(month(dep_time) == 1) %>%
mutate(dep_hour = update(dep_time, yday = 1)) %>%
ggplot(aes(dep_hour)) +
geom_freqpoly(binwidth = 300)
# 4月
flights_dt %>% filter(month(dep_time) == 4) %>%
mutate(dep_hour = update(dep_time, yday = 1)) %>%
ggplot(aes(dep_hour)) +
geom_freqpoly(binwidth = 300)
# 7月
flights_dt %>% filter(month(dep_time) == 7) %>%
mutate(dep_hour = update(dep_time, yday = 1)) %>%
ggplot(aes(dep_hour)) +
geom_freqpoly(binwidth = 300)
# 10月
flights_dt %>% filter(month(dep_time) == 10) %>%
mutate(dep_hour = update(dep_time, yday = 1)) %>%
ggplot(aes(dep_hour)) +
geom_freqpoly(binwidth = 300)
# あまり違いはなさそう

# 2.
flights_dt %>% mutate(sched_dep_time = update(sched_dep_time, yday = 1),
dep_time = update(dep_time, yday = 1)) %>%
ggplot() +
geom_point(aes(sched_dep_time, dep_time))

flights_dt %>% mutate(sched_dep_time = update(sched_dep_time, yday = 1),
dep_time = update(dep_time, yday = 1)) %>%
ggplot() +
geom_point(aes(sched_dep_time, dep_delay))
# sched_dep_timeの6:00のところで、数が増える
flights_dt %>% mutate(sched_dep_time = update(sched_dep_time, yday = 1),
dep_time = update(dep_time, yday = 1)) %>%
ggplot() +
geom_point(aes(dep_time, dep_delay))
# プロットがななめである

# 3.
flights_dt %>%
mutate(time_diff = arr_time – dep_time,
time_diff1 = ifelse(arr_time – dep_time < -1000,
arr_time – dep_time+1440, arr_time – dep_time )) %>%
ggplot() +
geom_point(aes(air_time, time_diff1))+
xlim(0,500)+ylim(0,500)
# 地上走行の時間があるので、air_timeより時間差の方が長い
# しかし、経度の差が大きいと、時差の影響で、時間差の方が短くなることがある

# 4.
flights_dt %>% group_by(hour(dep_time)) %>%
mutate(mean_dep_delay = mean(dep_delay)) %>%
ggplot() +
geom_line(aes(hour(dep_time), mean_dep_delay))

flights_dt %>% group_by(hour(sched_dep_time)) %>%
mutate(mean_dep_delay = mean(dep_delay)) %>%
ggplot() +
geom_line(aes(hour(sched_dep_time), mean_dep_delay))

# 一日の中では、時刻が遅くなるにつれて遅延も大きくなる
# sched_dep_Timeを使うべきである
# dep_timeの中にも遅延情報が入っており、一次独立でないからである

# 5.
flights_dt %>%
mutate(sched_dep_time = update(sched_dep_time, yday = wday(dep_time))) %>%
group_by(mday(sched_dep_time)) %>%
mutate(mean_dep_delay = mean(dep_delay)) %>%
ggplot() +
geom_line(aes(mday(sched_dep_time), mean_dep_delay))
# 土曜に出発すべき

# 6.
ggplot(diamonds) +
geom_freqpoly(aes(carat))

ggplot(flights_dt) +
geom_freqpoly(aes(sched_dep_time))

# 似ていないから何とも言えない

# 7.
flights_dt %>%
mutate(minute = minute(dep_time),
delay_flg = ifelse(dep_delay>0,1,0)) %>%
group_by(minute) %>%
select(minute, delay_flg) %>%
count() %>%
ggplot() +
geom_line(aes(minute, n))

13.4.4 まとめ

# 1.
# monthには30日と31日があり、定まらないため

# 2. Excelの式で説明?
# days(overnight * 1)
# := days(if(arr_time < dep_time, 1, 0) * 1)
# := days(if(arr_time < dep_time, 1 * 1, 0 * 1))
# := days(if(arr_time < dep_time, 1, 0))

# 3.
rep(ymd(“2015-1-1”),12) + months(c(0:11))

# 4.
cal_old <- function(birthday){
(ymd(birthday) %–% today()) %/% years(1)
}
# 中括弧のところでCtrl-Enterを押さないと
# object ‘birthday’ not foundのエラーが出る
cal_old(“19180401”)

# 5.
(today() %–% (today() + years(1)))/months(1)
# うまく働いているように見えるが。。。

15章 関数

15.2 関数を書くべきとき

# 1.
rescale01(TRUE)
range(TRUE)
# TRUE単体だとrangeが1~1になり、零割となる
# ただし、TRUEはほかの数と一緒になれば引数になりえる
rescale01(c(TRUE,FALSE))

# 欠損値を含むようにX1を作る
x1 <- c(1:10, NA, 12:20)
rescale01 <- function(x) {
rng <- range(x, na.rm = FALSE, finite = TRUE)
(x – rng[1]) / (rng[2] – rng[1])
}
rescale01(x1)
# 別に何も起こらないようにみえる

# 2.
rescale01 <- function(x) {
rng <- range(x, na.rm = TRUE)
ifelse((x – rng[1]) == (rng[2] – rng[1]),
1,
ifelse(-(x – rng[1]) == (rng[2] – rng[1]),
-1,
(x – rng[1]) / (rng[2] – rng[1])
)
)
}
rescale01(x)

# 3.
# mean(is.na(x))は、欠損値の比率を算出する
NARate <- function(x) {
mean(is.na(x))
}
y <- c(1:10, NA, 12:20)
NARate(y)

# x / sum(x, na.rm = TRUE)は、割合を示す
Rate <- function(x) {
x / sum(x, na.rm = TRUE)
}
Rate(y)

# sd(x, na.rm = TRUE) / mean(x, na.rm = TRUE)は
#標準偏差と平均の比であるが、何を意味するかは難しい
SDvsMEAN <- function(x) {
sd(x, na.rm = TRUE) / mean(x, na.rm = TRUE)
}
SDvsMEAN(y)

# 4.
variance <- function(x) {
var(x, na.rm = TRUE)
}
variance(y)

skewness <- function(x) {
n <- length(x)
v <- var(x)
m <- mean(x)
third.moment <- (1/(n – 2)) * sum((x – m)^3)
third.moment/(var(x) ^ (3/2))
}
skewness(rnorm(10))

# 5.
a <- c(1, NA, TRUE, NA)
b <- c(NA, 2, FALSE, NA)
sameNALoc <- function(x, y) {
if(length(x) != length(y)) {
0
}
else {
sum(is.na(x) & is.na(y))
}
}
sameNALoc(a, b)

# 6.
is_directory <- function(x) file.info(x)$isdir
is_readable <- function(x) file.access(x, 4) == 0

is_directory(“C:\\”)
is_readable(“C:\\”)
# ファイルシステムの状態を取得する
# 関数を作る際などに、あらかじめ状態を取得することで、エラーを回避できる

# 7.
# 問題の趣旨が分からないので、略

15.3 関数は人間とコンピュータのためのもの

#1.
f1 <- function(string, prefix) {
substr(string, 1, nchar(prefix)) == prefix
}
# 機能:接頭辞を確認する
# 名前:check_prefix

f2 <- function(x) {
if(length(x) <= 1) return(NULL)
x[-length(x)]
}
f2(c(1,2,3,4,5))
# 機能:ベクトルの最終要素を削除する
# 名前:rm_last_element

f3 <- function(x, y) {
rep(y, length.out = length(x))
}
# 機能:xベクトルのサイズ分だけyの要素を繰り返し出力する
# 名前:adjust_to_the_same_length
f3(c(1,2,3), c(4,5,6,7))
rep(c(4,5,6,7), length.out = 15)

# 2.
# 略

# 3.
?rnorm
?MASS::mvrnorm
rnorm(5, mean = 100)
MASS::mvrnorm(100, c(2,3), matrix(c(1,2,3,4), ncol = 2))
# rmormは1次元、mvrnormは多次元の乱数を発生する
# 一貫性の意味が分からない

# 4.
# norm_r(), norm_d()の方が接頭語が頭にあるので、自動補完が働く
# normとrやdは意味が分かれているので、アンダスコアを使って
# 分離している点も良い

15.4.3 コードのスタイル

# 1.
?’if’
?’ifelse’
# ifは条件文を括弧でくくる。ifelseは条件文の後はカンマである。
# ifはサブルーチンであり、文中で処理をして値を返さない。
# ifelseは関数であり、値を返す
# ifはNAでエラーになる、ifelseはNAを返す

# 2.
library(lubridate)
lubridate::now()
greetings <- function() {
h <- lubridate::hour(lubridate::now())
if (h < 10) {
“good morning”
} else if(h < 18) {
“good afternoon”
} else {
“good evening”
}
}
greetings()

# 3.
fizzbuzz <- function(x) {
s <- “”
if (x %% 3 == 0) {
s <- “fizz”
}
if (x %% 5 == 0) {
s <- stringr::str_c(s, “buzz”)
}
if (s == “”) x else s
}

fizzbuzz(15)

# 4.
cut(temp ,
breaks = c(-100, 0, 10, 20, 30, 100),
labels = c(“freezing”, “cold”, “cool”, “warm”, “hot”)
)

sense <- function (temp) {
cut(temp ,
breaks = c(-100, 0, 10, 20, 30, 100),
labels = c(“freezing”, “cold”, “cool”, “warm”, “hot”)
)
}

sense(10)

# 5.
switch(2,
`1` = 1 + 2,
`2` = 3 + 4,
`3` = 5 + 6
)
# “で囲えば使える

# 6.
switch(e,
a = ,
b = “ab”,
c = ,
d = “cd”
)
# Error: object ‘e’ not foundになる
?switch

15.5.4 遅延評価

# 1.
commas <- function(…) stringr::str_c(…, collapse = “, “)
letters <- letters[1:10]
commas(letters, collapse = “-“)
# collapseが…でcommasに渡され、
# 内部でstr_cにcollapseが複数わたされることになるため

# 2.
rule(“Title”, pad = “-+”)
# 2文字以上だと次の行に達してしまうため
# 下行のようにwidthを変更する
rule <- function(…, pad = “-“) {
title <- paste0(…)
width <- (getOption(“width”) – nchar(title)) %/% nchar(pad) – 5
cat(title, ” “, stringr::str_dup(pad,width), “\n”, sep = “”)
}

# 3.
?mean
# 極端に小さかったり大きかったりする数値を除外して平均をとる
# 運動競技などの採点で、審査員のばらつきの影響を少なくする

# 4.
cor
?cor
# corは相関係数を計算する
# 正規分布している場合はピアソン、離散値等の場合には
# ケンドール又はスピアマンの順位相関係数を使う
# デフォルトはピアソンである

16章 ベクトル

16.3.4 欠損値

# 1.
is.finite(NA)
!is.infinite(NA)
# NAを評価した時、is.finiteはFALSE、!is.infiniteはTRUEになる

# 2.
dplyr::near
# function (x, y, tol = .Machine$double.eps^0.5)
# {
#   abs(x – y) < tol
# }

# 3.
# 略

# 4.
# floor 切り捨て
# ceiling 切り上げ
# trunc 0に近い方へ切り捨て
# round IEEEの丸め

# 5.
search()
library(readr)
help(package = readr)
# 論理ベクトル:parse_logical
# 整数ベクトル:parse_integer
# 実数ベクトル:parse_number

16.4.5 要素抽出

# 1.
# mean(is.na(x)) は、ベクトルXの中でNAが出現する割合を返す
sum(!is.finite(x))
# sum(!is.finite(x))は、NAの数を数えてはくれるが無限大との区別はつかない

# 2.
?is.vector()
# ベクトルならTRUEを返す
?is.atomic
is.atomic(NULL)
# アトミックベクトルでないNULLに対し、is.atomic(NULL)でTRUEをかえすため

# 3.
?setNames
# Usage: setNames(object = nm, nm)
?purrr::set_names
# Usage: set_names(x, nm = x, …)
# 前者は引数が2つ、後者は2つ以上である

# 4.
x <- c(-4:0, NA, 0:4)

# a.
# [と[[のどちらでもよいが、ベクトルの要素にアクセスする場合は[が一般的か
x[length(x)]

# b.
x[1:(length(x)%/%2)*2]

# c.
x[-length(x)]

# d.
x[x %% 2 == 0 & !is.na(x)]

# 5.
x <- c(-2:5)
x
x[-which(x > 0)]
x[x <= 0]
# whichはリストから除去する要素のベクトルを与える
# x<=0はTRUE, FALSEのベクトルを与える
# 結果は 同じに見えるけど…

# 6.
x[10]
# NAになる
x<- c(abc = 1, def = 2, xyz = 5)
x[“xxx”]
# NAになる

16.5.3 調味料のリスト

# 1.

# 2.
library(tidyverse)
df <- tibble(
x = runif(5),
y = rnorm(5)
)
str(df)
str(df[1])
# 主な相違点はtibbleには列に名前がある

16.7.3 tibble

# 1.
hms::hms(3600)
# 01:00:00
x <- hms::hms(3600)

str(x)
# Classes ‘hms’, ‘difftime’  atomic [1:1] 3600
# ..- attr(*, “units”)= chr “secs”

typeof(x)
# [1] “double”

attributes(x)
# $units
# [1] “secs”
#
# $class
# [1] “hms”      “difftime”

# 2.
tb <- tibble::tibble(x = 1:5, y = c(1, 2))
# Error: Column `y` must be length 1 or 5, not 2
# エラーが生じて作れない

# 3.
tb <- tibble::tibble(x = list(1, 2, 3))
tb
# できる

17章 purrrでイテレーション

17.2 forループ

# 1.
# a.
mtcars
output <- vector(“double”, ncol(mtcars))
for (i in seq_along(mtcars)) {
output[[i]] <- mean(mtcars[[i]])
}
output

# b.
nycflights13::flights
output <- vector(“character”, ncol(nycflights13::flights))
for (i in seq_along(nycflights13::flights)) {
output[[i]] <- typeof(nycflights13::flights[[i]])
}
output

# c.
iris
output <- vector(“integer”, ncol(iris))
for (i in seq_along(iris)) {
output[[i]] <- length(unique(iris[[1]]))
}
output

# d.
m <- c(-10, 0, 10, 100)
output <- array(0, dim = c(10, length(m)))
for(i in 1:dim(output)[2]) {
for(j in 1:dim(output)[1]) {
output[[(i-1)*(dim(output)[[1]]) + j]] <- rnorm(1, mean = m[i])
}
}
output

# 2.
out <- “”
for (x in letters) {
out <- stringr::str_c(out, x)
}

paste(letters, collapse = “”)

#—

x <- sample(100)
sd <- 0
for (i in seq_along(x)) {
sd <- sd + (x[i] – mean(x)) ^ 2
}
sd <- sqrt(sd / (length(x) – 1))
sd

sd(x)

#—

x <- runif(100)
out <- vector(“numeric”, length(x))
out[1] <- x[1]
for (i in 2:length(x)) {
out[i] <- out[i-1] + x[i]
}
out

cumsum(x)

# 3.
# 略

# 4.
x <- runif(100000)
output <- vector(“integer”, 0)
for (i in seq_along(x)) {
output <- c(output, lengths(x[[i]]))
}
output

output <- vector(“integer”, length(x))
for (i in seq_along(x)) {
output[[i]] = length(x[[i]])
}
output

# これくらいだ10秒と1秒くらいの差になる

17.3.4 シーケンス長不明

# 1.
files <- dir(“C:\\Temp\\CSV\\”, pattern = “\\.csv$”, full.names = TRUE)
out <- vector(“list”, length(files))
out1 <- vector(“list”, 1)
for (i in seq_along(files)) {
out[[i]] <- read.csv(files[[i]])
out1 = merge(out1, out[[i]])
}
out
out1

# 2.
x <- rnorm(5)
output <- vector(“integer”, length(x))
for (nm in names(x)) {
output[[nm]] = 1
}
output
# xに名前がないとループされない

x <- rnorm(5)
output <- vector(“integer”, length(x))
output
names(x) <- c(“one”,”two”, “tree”)
x
for (nm in names(x)) {
output[[nm]] = 1
}
output
names(x)
# xに名前が無いと<NA>として処理される
# 名前が無いところが複数でも1つの<NA>となる

x <- rnorm(5)
output <- vector(“integer”, length(x))
output
names(x) <- c(“one”,”two”, “three”, “one”, “two”)
x
for (nm in names(x)) {
output[[nm]] = 1
}
output
names(x)
# 名前に重複がある場合、重複分は削除されて1回しか処理されない

# 3.
output <- vector(“double”, length(iris))
names(output) <- names(iris)
show_mean <- for (i in seq_along(iris)) {
if (typeof(iris[[i]]) == “double”) {
output[[i]] = mean(iris[[i]])
}
}
# 縦に出力する方法はちょっとわからない

# 4.
trans <- list(
disp = function(x) x * 0.0163871,
am = function(x) {
factor(x, labels = c(“auto”, “manual”))
}
)
for (var in names(trans)) {
mtcars[[var]] <- trans[[var]](mtcars[[var]])
}
# dispは、排気量をcu inからl(リットル)に変換する
# amは、0、1をそれぞれオートマとマニュアルに置き換える

17.4 forループと関数型

# 1.
?apply
# 第2の例:apply(x, 2, sort)
# sortなので
# for (i in 1:(length(x)-1)) {
#   for (j in 2:length(x)) {
#   }
# }
# ということを言わせたいのか?

# 2.

df <- tibble(
a = rnorm(10),
b = rnorm(10),
c = c(“a”, “b”, “c”, “d”, “e”, “f”, “g”, “h”, “i”, “j”),
d = rnorm(10)
)

col_summary <- function(df, fun) {
output <- vector(“double”, length(df))
for (i in seq_along(df)) {
if(is_numeric(df[[i]]) == TRUE) {
output[i] <- fun(df[[i]])
} else {
output[i] <- 0
}
}
output
}
col_summary(df, mean)

17.5.2 基本R

# 1.
#  a.
map_dbl(mtcars, mean)

#  b.
nycflights13::flights
map_chr(nycflights13::flights, typeof)

#  c.
iris
map_int(map(iris, unique), length)

# d.
m <- c(-10, 0, 10, 100)
map(m, rnorm, n = 10)

# 2.
as.character(map(iris, attr,”class”))
# str(iris)でもSpecisのみfactorであることが分かるが、
# 単一ベクトル化の方法は分からない

# 3.
map(1:5, runif)
# 要素のそれぞれにマッピングする
# 1:5が1~5の要素を作り、runifは、その数値分の乱数を返すので、
# 都合15個の乱数ができる

# 4.
map(-2:2, rnorm, n = 5)
# -2~2までのそれぞれの平均を持つ5個ずつの乱数が返る

map_dbl(-2:2, rnorm, n = 5)
# 動かない。
# 結果がベクトルにならないからである

# 5.
x <- list(list(1, 2, 3), list(4, 5, 6), list(7, 8, 9))
map(x, function(x) lm(x[[1]] ~ x[[2]], data = x))
# 意味がよくわからんが、こんな感じか?
# mpg, wtはmtcarsというData Frameにしかないので、そこを直す
# lmは線形近似する関数

17.9.2 reduceとaccumulate

# 1.
my_every <- function (.x, .p, …)
{
.p <- as_mapper(.p, …)
for (i in seq_along(.x)) {
val <- .p(.x[[i]], …)
if (rlang::is_false(val))
return(FALSE)
if (anyNA(val))
return(NA)
}
TRUE
}
x <- list(1:5, letters, list(10))
x %>% my_every(is_vector)

# 2.

df <- tibble(
a = rnorm(10),
b = rnorm(10),
c = rnorm(10),
d = rnorm(10)
)
col_summary <- function(df, fun) {
output <- vector(“double”, length(df))
for (i in seq_along(df)) {
if (typeof(df[[i]]) == “double”) {
output[i] <- fun(df[[i]])
} else {
output[i] <- NA
}

}
output
}
col_summary(mtcars, mean)
mtcars
typeof(iris[[1]])

# 3.
col_sum3 <- function(df, f) {
is_num <- sapply(df, is.numeric)
df_num <- df[, is_num]
sapply(df_num , f)
}
df <- tibble(
x = 1:3,
y = 3:1,
z = c(“a”, “b”, “c”)
)
col_sum3(df, mean)
col_sum3(df[1:2], mean)
col_sum3(df[1], mean)
col_sum3(df[0], mean)
df[0]
sapply(df[0], is.numeric)
# df[0]の戻り値が”# A tibble: 3 x 0″であり、
# sapply(df[0], is.numeric)が”named list()”と評価され、
# TRUE、FALSEが戻り値でないこと

18章 modelrを使ったモデルの基本

18.2 単純なモデル

# 1.
sim1a <- tibble(
x = rep(1:10, each = 3),
y = x * 1.5 + 6 + rt(length(x), df = 2)
)

sim1a_mod <- lm(y ~ x, data = sim1a)
coef(sim1a_mod)
sim1a_mod$coefficients[1]

ggplot(sim1a, aes(x, y)) +
geom_point(size = 2, color = “gray30”) +
geom_abline(intercept = sim1a_mod$coefficients[1],
slope = sim1a_mod$coefficients[2])
# 大きく離れた異常値があると、それに引きずられると言わせたいのかな?

# 2.
measure_distance <- function(mod, data) {
diff <- data$y – make_prediction(mod, data)
mean(abs(diff))
}

sim1a <- tibble(
x = rep(1:10, each = 3),
y = x * 1.5 + 6 + rt(length(x), df = 2)
)

best <- optim(c(0, 0), measure_distance, data = sim1a)

# Error in make_prediction(mod, data) :
# could not find function “make_prediction”
# となってしまいわからない

# 3.
model1 <- function(a, data) {
a[1] + data$x * a[2] + a[3]
}
# a[1]とa[3]が一次独立でないので、解が不定になる

18.3.2 残差

# 1.
sim1a_mod <- loess(y ~ x, data = sim1)
sim1a_mod

grid_a <- sim1 %>%
data_grid(x)

grid_a <- grid_a %>%
add_predictions(sim1a_mod)
grid_a

ggplot(sim1, aes(x, resid)) +
geom_point(aes(y = y)) +
geom_line(
aes(y = pred),
data = grid_a,
colour = “red”,
size = 1
)

ggplot(sim1, aes(x, y)) +
geom_point(aes(y = y)) +
geom_smooth(aes(y = y))

# geom_smoothもloessを使っているのでほぼ同じ

# 2.
# データフレームadd_predictionは、
# 入力データに単一の新しい列.predを追加します。
# spread_predictionsはモデルごとに1つの列を追加します。
# gather_prectionsは、2つの列.modelと.predを追加し、
# 各モデルに対して入力行を繰り返します。

# 3.
?geom_ref_line

# 参照線を追加します(ggplot2)。
ggplot(sim1, aes(x, resid)) +
geom_ref_line(h = 0) +
geom_point() +
geom_ref_line(h = -2, size=3, colour = “red”) +
geom_ref_line(h = 2, size=3, colour = “red”)

# ある閾値をこえる残差が分かるから?

# 4.

# 予測が妥当なら通常は正規分布になるはずだからである
# 利点:モデルの如何にかかわらず、統一的に評価できる
# 欠点:定義域に対する分布の状況が隠されてしまう

18.4.4 変換

# 1.
sim2
mod2 <- lm(y ~ x – 1, data = sim2)

grid <- sim2 %>%
data_grid(x) %>%
add_predictions(mod2)
ggplot(sim2, aes(x)) +
geom_point(aes(y = y)) +
geom_point(
data = grid,
aes(y = pred),
color = “red”,
size = 4
)
sim2
mod2
# 別に何も起こらないように見える
# 繰り返しの方法が分からない

# 2.
model_matrix(sim3, y ~ x1 + x2)
model_matrix(sim3, y ~ x1 * x2)

model_matrix(sim4, y ~ x1 + x2)
model_matrix(sim4, y ~ x1 * x2)

# * の方が組合せが多いから

# 3.
mod1 <- lm(y ~ x1 + x2, data = sim3)
fnc1 <- function(.x) {
predict(mod1, .x)
}
data <- tribble(
~ x1, ~ x2,
1, “a”,
1, “b”,
2, “a”,
2, “b”
)
fnc1(data)

mod2 <- lm(y ~ x1 * x2, data = sim3)
fnc2 <- function(.x) {
predict(mod2, .x)
}
fnc2(data)

# 4.
sim4
sim41 <- sim4 %>%
add_residuals(mod1)
ggplot(sim41, aes(resid)) +
geom_freqpoly(binwidth = 0.5)

sim42 <- sim4 %>%
add_residuals(mod2)
ggplot(sim42, aes(resid)) +
geom_freqpoly(binwidth = 0.5)
# これで良いのかわからない。うまく書けない

19章 モデル構築

19.2.2 さらに複雑なモデル

# 1.
# 明るい縦の縞は、カウント数が多いことを示している
# lcaratの区切りの良い数(-1, 0, 1)などに有ることから、
# ダイヤモンドは、区切りの良いカラット数でカットされることが多い

# 2.
# log(price) = a_0 + a_1 * log(carat)
# log(price) – log(exp(a_0)) = a_1 * log(carat)
# log(price/exp(a_0)) = log(carat ^ a_1)
# price/exp(a_0) = carat ^ a_1
# price = exp(a_0) * (carat ^ a_1)
# つまり、価格をカラットの指数関数でモデルすることになる

# 3.
diamonds2 %>%
filter(lresid2 == max(diamonds2$lresid2))

# caratが0.340、cutはFair、clarityはI1でかなり粗悪なダイヤモンドであり、
# 値段の付け間違えか、詐欺用途的?

diamonds2 %>%
filter(lresid2 == min(diamonds2$lresid2))

# caratが2.46、cutがPremium、clarityがSI2で、大きいものの、普通の出来
# 大きいダイヤが欲しければ、買っても良いことになる

# 4.
# 簡易的には、判断根拠はP.340のlcaratとlresid2のグラフになる
# このグラフは、モデリングに用いたデータを使って作成しているので、
# 本来なら評価用のデータをあらかじめ別に確保しておかなければならない
#
# 全体的には、ほぼ、平らに一様分布している
# lcaratの小さい部分では正側に大きい残差があり、
# 大きい部分では負側に大きい残差がある
# lcaratの大きい部分で正側に残差が出ない理由の一つは、
# $19,000ドルを超えないデータを用いていることが一因である
# よって、lcaratの高い領域では、判断が付きかねるので、
# 低い領域、具体的には1カラット以下のダイヤモンドを
# 買うに当たっては、信用しても良いのではないかと考える

19.3.4 年間の日:別のアプローチ

# 1.
日曜日ということ以外よく分からない

# 2.
daily %>%
add_residuals(mod1, “resid”) %>%
top_n(3, resid)
# クリスマスで人の移動が多いのかな?
# 一般化するという意味が分からない

# 3.
#https://jrnold.github.io/r4ds-exercise-solutions/model-building.htmlを見た
daily <- daily %>%
mutate(
wday2 =
case_when(
wday == “土” & term == “summer” ~ “Sat-summer”,
wday == “土” & term == “fall” ~ “Sat-fall”,
wday == “土” & term == “spring” ~ “Sat-spring”,
TRUE ~ as.character(wday)
)
)

mod4 <- lm(n ~ wday2, data = daily)
daily %>%
gather_residuals(combination = mod2, with_term = mod4) %>%
ggplot(aes(date, resid, color = model)) +
geom_line(alpha = 0.75)
# あらゆる組合せをしたモデルの方が良好である。特に夏期。

# 4.
# 祝日が分からないので略

# 5.
daily <- daily %>%
mutate(
month = month(date)
)
mod5 <- lm(n ~ wday * month, data = daily)
daily %>%
gather_residuals(with_term = mod2, whth_month = mod5) %>%
ggplot(aes(date, resid, color = model)) +
geom_line(alpha = 0.75)
# termの変わり目が月の途中(6/5、8/25)にあるため

# 6.
mod6 <- lm(n ~ wday + ns(date,5), data = daily)
daily %>%
gather_residuals(with_term = mod2, whth_ns = mod6) %>%
ggplot(aes(date, resid, color = model)) +
geom_line(alpha = 0.75)
# nsはスプライン補間であり、祝日などの特異日や
# 自由度数の合わない場合には精度が悪い

# 7.
library(lubridate)
daily2 <- flights %>%
mutate(
date = make_date(year, month, day),
wday = wday(date, label = TRUE))

daily2 <- daily2 %>%
filter(wday == “日”) %>%
mutate(
dep_time1 = dep_time %/% 100,
distance1 = distance %/% 100 * 100) %>%
select(dep_time = dep_time1, distance = distance1) %>%
group_by(dep_time, distance)

ggplot(daily2, aes(distance, dep_time)) +
geom_hex(bins = 20)
# 遠くだからといって、夜の便が多くなるということは認められない

# 8.
# ファクタの水準を設定する関数
wdayMonStart <- function(wday) {
fct_relevel(wday,”月”,”火”,”水”,”木”,”金”,”土”)
}
# テスト
daily <- daily %>%
mutate(wday = wday(date, label = TRUE))
ggplot(daily, aes(wdayMonStart(wday), n)) +
geom_boxplot()
fct_relevel

20章 purrrとbroomによる多数のモデル

20.2.4 モデル品質

# 1.
gapminder <- gapminder %>%
mutate(newyear = year – mean(year))
by_country_new <- gapminder %>%
group_by(country, continent) %>%
nest()
country_model_new <- function(df) {
lm(lifeExp ~ poly(newyear, 2), data = df)
}
models_new <- map(by_country_new$data, country_model_new)
by_country_new <- by_country_new %>%
mutate(model_new = map(data, country_model_new))
by_country_new <- by_country_new %>%
mutate(
resids = map2(data, model_new, add_residuals)
)
resids_new <- unnest(by_country_new, resids)
resids_new

resids_new %>%
ggplot(aes(year, resid)) +
geom_line(aes(group = country), alpha = 1 / 3) +
geom_smooth(se = FALSE)

resids_new %>%
ggplot(aes(year, resid, group = country)) +
geom_line(alpha = 1 / 3) +
facet_wrap(~continent)
# Americas, Asiaでは上に凸の形状が少なくなった
# しかし、Africaが合わないのは相変わらずである

# 2.
library(ggbeeswarm)
glance %>%
ggplot(aes(continent, r.squared)) +
geom_quasirandom()
# 下に凸の正規分布的なグラフが書けた

# 3.
glance <- by_country %>%
mutate(glance = map(model, broom::glance)) %>%
unnest(glance)

# 分かりません

20.4.4 名前付きリストからリスト列を作る

# 1.
# 下記等を実行して、出てきた関数一覧から探すくらいしか思いつかない
ls(“package:tidyr”)

# 2.
# 独学なので、ブレインストーミングできない

# 3.
mtcars %>%
group_by(cyl) %>%
summarize(q = list(quantile(mpg))) %>%
unnest()

mtcars %>%
ggplot(aes(mpg)) +
geom_histogram(binwidth = 1) +
facet_wrap(~cyl)

mpg6 <- mtcars %>%
filter(cyl == 6) %>%
select(mpg)
mpg6[order(mpg6$mpg),]

# 例えば、quantileでは、cyl == 6に対するmpgの第1四分位が18.6と
# 計算されるが、mtcarsの中に、(cyl, mpg) = (6, 18.6)というサンプルは無い

# 4.
mtcars %>%
group_by(cyl) %>%
summarize_each(funs(list))

?funs

# funsは関数呼び出しのリストを作る。
# 他の関数に入力するための関数の名前付きリストを生成するための
# 柔軟な方法を提供するので役立つ

20.5.2 入れ子解消

# 1.
# 本文P.371にも書いてあるように、欠損値等があった場合に
# filterで取り除くことができるから

# 2.
number <- c(1, 2, 3)
id <- c(10, 20, 30)
name <- c(“次郎左衛門”, “平兵衛”, “次郎吉”)
x <- data.frame(Number = number, ID = id, Name = name)
x
sapply(x, class)
# 一般的な型はnumericかな?
# 一次元と多次元の違い?

21章 Rマークダウン

21.2 Rマークダウン

# 1.


title: “My Mame”
output: html_document

# 学歴
* __2001.04.01-2004.03.31__ 学位1
* __2004.04.01-2007.03.31__ 学位2
* __2007.04.01-2010.03.31__ 学位3

# 職歴
1. ___2010.04.01-2013.03.31___ 職種1
1. ___2013.04.01-2016.03.31___ 職種2
1. ___2016.04.01-2019.03.31___ 職種3

*   Bulleted list item 1

*   Item 2

* Item 2a

* Item 2b

# 2.
## a. 脚注追加
脚注 [^1]
[^1]: これが脚注です。脚注番号はうまく働かないみたい

## b. 水平区切り線追加
***

## c. 引用追加
> 引用追加

## 3.

“`{r}

“`

“`{r}

“`

“`{r by-name}
mtcars[1:5, 1:10]

knitr::kable(
mtcars[1:5, ],
caption = “A knitr kable.”
)
“`

21.4.6 インラインコード

## 1.
“`{r}
knitr::opts_chunk$set(
echo = FALSE
)
library(ggplot2)
“`

ダイヤモンドの大きさは、カットに対し、次のように分布します、
“`{r}
ggplot(diamonds, aes(carat)) +
geom_histogram(binwidth = 0.01) +
facet_wrap(~ cut)
“`

また、色と透明度では次のように分布します。
“`{r}
ggplot(diamonds, aes(carat)) +
geom_histogram(binwidth = 0.01) +
facet_wrap(~ color)
ggplot(diamonds, aes(carat)) +
geom_histogram(binwidth = 0.01) +
facet_wrap(~ clarity)
“`

# 3. 略

# 4. 略

22章 ggplot2でコミュニケーションのためのグラフ作成

22.2 ラベル

# 1.
ggplot(mpg, aes(displ, hwy)) +
geom_point(aes(color = class)) +
geom_smooth(se= FALSE) +
labs(
title = “燃費と排気量の関係”,
subtitle = “ただし、燃費はガロンとマイルです”,
caption = “データはfueleconomy.govから”,
x = quote(lim(f(x), x %->% 0)),
y = quote(integral(f(x)*dx, a, b))
)

# 2.
mpg[order(mpg$displ, decreasing = TRUE),]
# 車種の列はclassであることがわかる
library(modelr)
library(splines)
mod <- MASS::rlm(hwy ~ ns(displ, 2) * class, data = mpg)
mpg %>%
data_grid(displ = seq_range(displ, n = 10), class, hwy) %>%
add_predictions(mod) %>%
ggplot(aes(displ, pred, color = class)) +
geom_line() +
geom_point()
# 2seaterの低排気量域とcompctの大排気量域がうまくいかない

# 3.
mpg %>%
data_grid(displ = seq_range(displ, n = 10), class, hwy) %>%
add_predictions(mod) %>%
ggplot(aes(displ, pred, color = class)) +
geom_line() +
geom_point() +
labs(
title = “燃費と排気量のモデル”,
subtitle = “2seaterの低排気量域とcompctの大排気量域がうまくいっていません”,
caption = “データはfueleconomy.govから”,
x = “排気量[l]”,
y = “高速走行時の燃費[gallon/miles]”
)

22.3 アノテーション

# 1.
label <- tribble(
~ displ, ~ hwy, ~ label,
Inf,   Inf,    “++”,
Inf,  -Inf,    “+-“,
-Inf,   Inf,    “-+”,
-Inf,  -Inf,    “–”
)

ggplot(mpg, aes(displ, hwy)) +
geom_point() +
geom_text(aes(label = label),
data = label,
vjust=”top”,
hjust = “right”) +
geom_text(aes(label = label),
data = label,
vjust=”bottom”,
hjust = “right”) +
geom_text(aes(label = label),
data = label,
vjust=”top”,
hjust = “left”) +
geom_text(aes(label = label),
data = label,
vjust=”bottom”,
hjust = “left”)
# 2.
?annotate
ggplot(mpg, aes(displ, hwy)) +
geom_point() +
annotate(“text”, x = 4, y = 30, label = “Some text”)

# 3.
ggplot(mpg, aes(displ, hwy)) +
geom_point(aes(color = class)) +
geom_text(aes(label = model), data = best_in_class) +
facet_wrap(~ class)
# facetにしただけで、各ファセットに異なるラベルが付いた

# 4.
?geom_label
# label.padding: ラベルの周囲の余白の量
# label.r      : 箱の四隅の半径
# label.size   : 境界線の太さ(mm)

# 5.
?arrow()
# angle : 矢印の角度(degree)
# length: 長さ
# ends  : “last”, “first” 又は”both”により、矢印を付ける位置を定める
# type  : “open”又は”closed”。closedの場合は矢印の先頭を塗りつぶした三角形にする

22.4.3 スケールの入れ替え

# 1.
ggplot(df, aes(x, y)) +
geom_hex() +
scale_color_gradient(low = “white”, high = “red”, space = “Lab”) +
coord_fixed()

# 分からない。geom_hexが対応していないから?
# 色を付ける対象が明記されていないから?

# 2.
? scale_color_gradient
# 「名前、制限、分割、ラベルなどを制御するために
#  continuous_scaleに渡されるその他の引数」である
#
?labs
# 「新しい名前と値のペアのリスト」である

# 3.
library(modelr)
library(lubridate)

presidential %>%
mutate(id = 33 + row_number()) %>%
ggplot(aes(start, id)) +
geom_ref_line(v = ymd(19520101)) +
geom_ref_line(v = ymd(19560101)) +
geom_ref_line(v = ymd(19600101)) +
geom_ref_line(v = ymd(19640101)) +
geom_ref_line(v = ymd(19680101)) +
geom_ref_line(v = ymd(19720101)) +
geom_ref_line(v = ymd(19760101)) +
geom_ref_line(v = ymd(19800101)) +
geom_ref_line(v = ymd(19840101)) +
geom_ref_line(v = ymd(19880101)) +
geom_ref_line(v = ymd(19920101)) +
geom_ref_line(v = ymd(19960101)) +
geom_ref_line(v = ymd(20000101)) +
geom_ref_line(v = ymd(20040101)) +
geom_ref_line(v = ymd(20080101)) +
geom_ref_line(v = ymd(20120101)) +
geom_ref_line(v = ymd(20160101)) +
geom_point() +
geom_segment(aes(xend = end, yend = id)) +
scale_x_date(
NULL,
breaks = presidential$start,
date_labels = “‘%y”
) +
labs(
y = “第n代大統領”
) +
ggrepel::geom_label_repel(aes(label = name), alpha = 0.6) +
ggtitle(“第n代大統領と任期の関係”)

# 4.
ggplot(diamonds, aes(carat, price)) +
geom_point(aes(color = cut), alpha = 1/20) +
guides(color = guide_legend (
override.aes = list(alpha = 1, size = 4)))

出典

『R for Data Science』(Hadley Wickham、Garrett Grolemund 著、O’Reilly、Copyright 2017 O’Reilly Media、ISBN978-1-491-91039-9、日本語版『Rではじめるデータサイエンス』オライリー・ジャパン、ISBN978-4-87311-814-7

シェアする

  • このエントリーをはてなブックマークに追加

フォローする