5章 探索的データ分析

5.0 ライブラリーの読み込み

library("tidyverse")
library("viridis")
library("nycflights13")
library("lvplot")
library("ggbeeswarm")
library("ggstance")
library("gridExtra")

5.1 はじめに

練習問題はありません

5.2 質問

練習問題はありません

5.3 変動

練習問題1 : diamondsのx、y、zの分布を探索しなさい。また、関係性を考え直さい。

まずは、summary()を使って、基本的な要約統計量を確認します。yとzには外れ値があるように見えます。また、ダイアモンドの形状を表すx、yとzであるにも関わらず、いくつかのダイヤモンドはゼロの値を持っているようです。

diamonds %>% 
   select(x, y, z) %>% 
   summary()
   
       x                y                z         
 Min.   : 0.000   Min.   : 0.000   Min.   : 0.000  
 1st Qu.: 4.710   1st Qu.: 4.720   1st Qu.: 2.910  
 Median : 5.700   Median : 5.710   Median : 3.530  
 Mean   : 5.731   Mean   : 5.735   Mean   : 3.539  
 3rd Qu.: 6.540   3rd Qu.: 6.540   3rd Qu.: 4.040  
 Max.   :10.740   Max.   :58.900   Max.   :31.800  

視覚化することでよりわかりやすくなります。

x <- list("x","y","z")
plots <- x %>% 
  map( ~ ggplot(data = diamonds) + 
         geom_histogram(aes_string(x = .x), binwidth = 0.05) + ylab(.x)) 
grid.arrange(plots[[1]], plots[[2]], plots[[3]],
             ncol = 3, nrow = 1)

# これでも同じなので好きなのをどうぞ
# x <- c("x", "y", "z")
# plots <- list()
# for (i in seq_along(x)) {
#   plots[[i]] <-
#     ggplot(diamonds, aes_string(x = x[[i]])) +
#     geom_histogram(binwidth = 0.05) +
#     ylab(x[[i]])
# }
# 
# grid.arrange(plots[[1]], plots[[2]], plots[[3]],
#              ncol = 3, nrow = 1)

これだと分布の形状が確認しづらいので、四分位範囲が収まる範囲でxlim(0,10)を設定して視覚化してみるとだいぶ見通しが良くなります。これを見る限り、ある特定の閾値まで行くと、意図的に数値が丸められているようです。

x <- list("x","y","z")
plots <- x %>% 
  map( ~ ggplot(data = diamonds) + 
         geom_histogram(aes_string(x = .x), binwidth = 0.05) + ylab(.x) + xlim(0,10)) 
grid.arrange(plots[[1]], plots[[2]], plots[[3]],
             ncol = 3, nrow = 1)

外れ値も可視化しておきます。この図を見る限り、xは3~10、yも3~10、zも3~10くらいの範囲があれば、大部分は可視化して傾向を把握できそうです。

x <- list("y","z")
plots <- x %>% 
  map( ~ ggplot(data = diamonds, aes(x = x)) + 
         geom_jitter(aes_string(y = .x)) + ylab(.x))
grid.arrange(plots[[1]], plots[[2]],
             ncol = 2, nrow = 1)

3から10の範囲に絞って可視化します。ある閾値まで行くと、意図的に数値がスパイクする様子が先程と同様に確認できます。

diamonds %>% 
  select(x, y, z) %>% 
  gather(key = shape, val = val ,.) %>% 
  filter(val > 0, val < 10) %>% 
  ggplot(., aes(val, fill = shape)) + 
  geom_histogram(position = "identity", binwidth = 0.01) +
  scale_x_continuous(breaks = 1:10) + 
  facet_grid(. ~ shape, scales = "free_x")

オーバーラップさせると、xとyは似たような変動をしていることがわかります。おそらく相関も高いことが予想されます。

diamonds %>% 
  select(x, y, z) %>% 
  gather(key = shape, val = val ,.) %>% 
  filter(val > 0, val < 10) %>% 
  ggplot(., aes(val, fill = shape)) + 
  geom_histogram(position = "identity", binwidth = 0.01, alpha = 1/2) +
  scale_x_continuous(breaks = 1:10)

練習問題2 : 価格の分布を調べなさい。何らかの発見を見つけなさい。

summary()で要約量統計量を見てみると、平均が3933で、最大値は18823にもなるようです。

diamonds %>%
  select(price) %>% 
  summary()

     price      
 Min.   :  326  
 1st Qu.:  950  
 Median : 2401  
 Mean   : 3933  
 3rd Qu.: 5324  
 Max.   :18823 

分布を可視化してみると、価格が突然ジャンプしているように見えます。

diamonds %>% 
  select(price) %>% 
  ggplot(aes(x = price)) +
  geom_histogram(binwidth = 100) +
  scale_x_continuous(breaks = seq(0, 20000, by = 1000))

さらに細かく探索すると、1450~1550の間の値がないことがわかります。

diamonds %>% 
  select(price) %>% 
  filter(between(price, 1000, 2000)) %>% 
  ggplot(aes(x = price)) +
  geom_histogram(binwidth = 5) +
  scale_x_continuous(breaks = seq(1000, 2000, by = 50))

練習問題3 : 0.99カラットのダイヤモンドはいくつですか？1カラットはいくつですか？違いは何だと思いますか？

0.99カラットのダイアモンドは23個です。0.9~1.0にかけて、谷ができているように見えます。

diamonds %>%
  select(carat) %>% 
  filter(between(carat, 0.90, 1.00)) %>% # filter(carat >= 0.99, carat <= 1)
  count(carat)

# A tibble: 11 x 2
   carat     n
   <dbl> <int>
 1  0.9   1485
 2  0.91   570
 3  0.92   226
 4  0.93   142
 5  0.94    59
 6  0.95    65
 7  0.96   103
 8  0.97    59
 9  0.98    31
10  0.99    23
11  1     1558

可視化して確認してみます。caratはどうやら小数点が9に近づくに連れて、切り上げられてラベルされているように見えます。

diamonds %>%
  select(carat) %>% 
  filter(between(carat, 0, 1.00)) %>%
  count(carat) %>% 
  ggplot(aes(carat, n)) +
  geom_histogram(stat = "identity", binwidth = 5) + 
  scale_x_continuous(breaks = seq(0.2, 1.0, by = 0.05))

練習問題4 : coord_cartesian()とxlim()またはylim()を使った場合の違いを調べなさい。

coord_cartesian()は計算してからズームアップして範囲を切り取るイメージです。反対に、*lim()は、計算する前に除外して範囲を切り取るイメージです。したがって、coord_cartesian()は範囲外にも値がありますが、*lim()の方にはありません。

g1 <- ggplot(diamonds) +
  geom_histogram(mapping = aes(x = price), binwidth = 100) +
  coord_cartesian(xlim = c(1000, 5000), ylim = c(0, 3000))

g2 <- ggplot(diamonds) +
  geom_histogram(mapping = aes(x = price), binwidth = 100) +
  xlim(1000, 5000) + ylim(0, 3000)

grid.arrange(g1, g2, ncol = 2)

5.4 欠損値

練習問題1 : ヒストグラムの欠損値はどうなりますか？棒グラフの欠損値はどうなりますか？

警告メッセージが表示されます。17502件はNAなので、削除されました。値がわからないので、ビンに入れようがありません。

diamonds %>%
  select(carat) %>% 
  mutate(carat = if_else(carat > 1, NA_real_, carat)) %>% 
  ggplot(., aes(x = carat)) +
  geom_histogram()
  
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
 
警告メッセージ: 
Removed 17502 rows containing non-finite values (stat_bin). 

カテゴリカル変数の場合はNAとして、欠損値の数をカウントすることができます。

diamonds %>%
  mutate(color = if_else(runif(n()) < 0.3, NA_character_, as.character(color))) %>%
  ggplot() +
  geom_bar(mapping = aes(x = color))

練習問題2 : mean()とsum()のna.rm = TRUEの挙動を調べなさい。

これらの関数はNAがベクトルの要素として混じっていると、NAを返すようになっている。そのため、NAを除いて計算させるためには、na.rm = TRUEを付ける必要がある。

diamonds %>%
  select(price) %>% 
  mutate(price = if_else(runif(n()) < 0.3, NA_integer_, price)) %>%
  summarise(avg = mean(price),
            sum = sum(price))
            
# A tibble: 1 x 2
    avg   sum
  <dbl> <int>
1    NA    NA

diamonds %>%
  select(price) %>% 
  summarise(avg = mean(price, na.rm = TRUE),
            sum = sum(price, na.rm = TRUE))

# A tibble: 1 x 2
    avg       sum
  <dbl>     <int>
1 3930. 148320930

5.5 共変動

練習問題1 : キャンセルフライトとフライトの出発時刻の可視化を改善しなさい。

箱ひげ図も使うことができるが、ここはバイオリンプロットと併用してみます。

flights %>%
  select(dep_time, sched_dep_time) %>% 
  mutate(
    cancelled = is.na(dep_time),
    sched_hour = sched_dep_time %/% 100,
    sched_min = sched_dep_time %% 100,
    sched_dep_time = sched_hour + (sched_min / 60)
  ) %>%
  ggplot(aes(x = cancelled, y = sched_dep_time)) +
  geom_violin() + 
  geom_boxplot(width=.1, alpha = 1/2, fill = "#006E4F")

練習問題2 : ダイヤモンドの価格を予測する上で、どの変数が最も重要ですか。その変数はカットとどのように関連していますか？これら2つの関係を組み合わせると、低品質のダイヤモンドがより高価になるのはなぜですか？

まずは連続変数とpriceの関係を探索していく。散布図を見る限りpriceと関係がありそうなのは、caratとxであり、これらは値が大きなるにつれて、priceも大きくなる傾向がある。caratは直感的にもおかしくないが、xはダイヤモンドの寸法なので、x軸の大きさのみ、ダイヤモンドの価格と関係しているのかもしてない。

num_plt <- diamonds %>% 
  select(-price) %>% 
  keep(is.numeric) %>% 
  names() %>% 
  map(list()) 

plots <- num_plt %>% 
  map( ~ ggplot(data = diamonds, aes(y = price)) + 
         geom_point(aes_string(x = .x)) + xlab(.x))

grid.arrange(plots[[1]], plots[[2]], plots[[3]],plots[[4]],plots[[5]],plots[[6]],
             ncol = 3, nrow = 2)

次はカテゴリカル変数について見ていきます。colorはDよりもJのほうがpriceが高い傾向にあります。clarityはl1(最悪)よりもIF(最良)はpriceが低い傾向にあります。

fctr_plt <- diamonds %>% 
  select(-price) %>% 
  keep(is.factor) %>% 
  names() %>% 
  map(list()) 

plots <- fctr_plt %>% 
  map( ~ ggplot(data = diamonds, aes(y = price)) + 
         geom_boxplot(aes_string(x = .x)) + xlab(.x))

grid.arrange(plots[[1]], plots[[2]], plots[[3]],
             ncol = 3, nrow = 1)

cutとの組み合わせを確認します。cutとcaratには否定的な関係があります。最大のカラットダイヤモンドにはFair(最低)のカットがあります。大きいダイヤモンドほど低品質のカットでも販売できますが、小さいダイヤモンドほど良いカットが必要なのでしょうか。

diamonds %>% 
  ggplot(aes(cut, carat, col = price)) + 
  geom_jitter(alpha = 0.2)

練習問題3 : {ggstance}パッケージを利用し、水平箱ひげ図を作成しなさい。coord_flip()の違いは何ですか？

まずは、geom_boxplot()とcoord_flip()を使って箱ひげを書きます。xにclassを使用して、座標を変換しています。{ggstance}のgeom_boxploth()では、素直にyにclassを使用して作図できます

g1 <- ggplot(data = mpg) +
  geom_boxplot(mapping = aes(x = reorder(class, hwy, FUN = median), y = hwy)) +
  coord_flip()

g2 <- ggplot(data = mpg) +
  geom_boxploth(mapping = aes(y = reorder(class, hwy, FUN = median), x = hwy))

grid.arrange(g1, g2, ncol = 2)

練習問題4 : {lvplot}パッケージを利用し、geom_lv()でprice vs cutの分布を表示しなさい。

geom_lv()は、箱ひげを改良した作図法で、よりも多い分位数を使って作図します。

ggplot(diamonds, aes(x = cut, y = price)) +
  geom_lv()

練習問題5 : geom_violin()、geom_histogram()、geom_freqpoly()を比較し、各方法の長所と短所を調べなさい。

geom_freqpoly()は、priceが最もカウント(y = ..density..で密度も可能)されている場所を特定したり、他のカテゴリとのカウント数の違いを見分けるのは簡単です。geom_histogram()も、positionを工夫すれば、似たようなことは可能です。

geom_violin()は、分布の要約値(歪度、中心値、分散など)を視覚的に区別するのは簡単です。しかし、カウントの多さなどは比較できません。

g1 <- ggplot(data = diamonds, mapping = aes(x = price)) +
  geom_freqpoly(mapping = aes(color = cut), binwidth = 500) +
  guides(colour = FALSE) 

g2 <- ggplot(data = diamonds, mapping = aes(x = price)) +
  geom_histogram(mapping = aes(fill = cut), position = "dodge", alpha = 0.8, binwidth = 500) +
  guides(fill = FALSE) 

g3 <- ggplot(data = diamonds, mapping = aes(x = cut, y = price)) +
  geom_violin(mapping = aes(fill = cut)) +
  guides(fill = FALSE) 

grid.arrange(g1, g2, g3, ncol = 3)

練習問題6 : 小さいデータの場合geom_jitter()は、連続変数とカテゴリカル変数の関係を調べるために便利です。{ggbeeswarm}パッケージは似た多くの方法を提供しています。それらをリスト化し、何をするのか簡単に説明してください。

• geom_quasirandom()：ジッタープロットとバイオリンプロットを組み合わせたプロットを作成します。method = c("tukey", "frowney", "smiley")があります。

• geom_beeswarm() ：バイオリンプロットに似たプロットを作成しますが、点をオフセットします。

g1 <- ggplot(data = mpg) +
  geom_quasirandom(
    mapping = aes(
      x = reorder(class, hwy, FUN = median),
      y = hwy
    ),
    method = "tukey"
  )

g2 <- ggplot(data = mpg) +
  geom_beeswarm(mapping = aes(
    x = reorder(class, hwy, FUN = median),
    y = hwy
  ))

grid.arrange(g1, g2,ncol = 2)

練習問題7 : color内のcutの分布をより明確に示すために、データセットをどのようにスケール変更すればよいか。

cut内の分布を明確に示すために、ここでは全体に占めるcutとcolorの組み合わせ数の割合で色をつけてみます。color:G, cut:idealが多いことがわかります。

diamonds %>%
  count(color, cut) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(mapping = aes(x = color, y = cut)) +
  geom_tile(mapping = aes(fill = prop)) +
  scale_fill_gradient(low = "#FFFFFF", high = "#006E4F")

練習問題8 : geom_tile()と{dplyr}パッケージを一緒に使用し、目的地や年月によって平均飛行遅延がどのように変わるかを調べなさい。

比較する数が多くなると、数字での比較よりも、視覚的に表現するほうがわかりやすい。傾向として、夏は平均遅延時間が長くなる傾向があり、特定のキャリア(CAE)は一貫して平均遅延時間が長いことがわかります。

flights %>%
  group_by(month, dest) %>%
  summarise(dep_delay = mean(dep_delay, na.rm = TRUE)) %>%
  ggplot(aes(x = factor(month), y = dest, fill = dep_delay)) +
  geom_tile() +
  scale_fill_gradient(low = "#FFFFFF", high = "#006E4F") +
  labs(x = "Month", y = "Destination", fill = "Avg Dep Delay")

練習問題9 : さきほどの例ではaes(x = color, y = cut)なく、aes(x = cut, y = color)のほうが、よいのでしょうか。

通常、カテゴリが多い変数をy軸に置くことをお勧めします。

練習問題10 : 条件付きの分布を箱ひげ図で要約する代わりに、度数分布多角形を使用できる。cut_width()とcut_number()を使うときに何を考慮する必要があるか。

cut_number()もcut_width()も分類分けする関数です。これらを使用する場合、カットする範囲でプロットの印象が大きく変わるため、カット範囲には注意しなければいけない。

g1 <- ggplot(data = diamonds,　mapping = aes(x = price, color = cut_number(carat, 5))) +
  geom_freqpoly() +
  labs(x = "Price", y = "Count", color = "Carat")

g2 <- ggplot(data = diamonds,mapping = aes(x = price, color = cut_width(carat, 1, boundary = 0))) +
  geom_freqpoly() +
  labs(x = "Price", y = "Count", color = "Carat")

grid.arrange(g1, g2, ncol = 2)

練習問題11 : priceで区分けしたcaratの分布を視覚化しなさい。

さきほど学んだcut_width()で2000単位でpriceを区切って視覚化してみます。

ggplot(diamonds, aes(x = cut_width(price, 2000, boundary = 0), y = carat)) +
  geom_boxplot() + coord_flip() + xlab("Price")

練習問題12 : 大きいダイヤモンドと小さいダイヤモンドの価格分布を比較しなさい。

大きいダイアモンド(5,6]もそもそも数が少なく、分布を把握するのが困難である一方で、[0,1]の分布は外れ値のような価格が高いものと低いものが混じっている分布のようです。

ggplot(diamonds, aes(x = cut_width(carat, 1, boundary = 0), y = price)) +
  geom_boxplot() + xlab("carat")

priceを離散化して、ヒートマップで表現してみます。全体的な傾向として、ポジティブな相関が確認でき、caratが小さいのに、priceが高いものや、そもそもcaratが大きいものは全体に占める割合は小さいようです。

diamonds %>% 
  select(carat, price) %>% 
  mutate(carat = cut_width(carat, .5, boundary = 0),
         price = cut_width(price, 1000, boundary = 0)) %>% 
  complete(carat, price, fill = list(n = 0L)) %>%
  group_by(carat, price) %>% 
  summarise(n = n()) %>% 
ggplot(., aes(x = carat, y = price, fill = n)) +
  geom_tile() +
  scale_fill_gradient(low = "#FFFFFF", high = "#006E4F") 

練習問題13 : カット、カラット、および価格の組み合わせ分布を視覚化するために、学んだ2つの手法を組み合わせなさい。

散布図で一気に視覚化しつつ、アルファブレンディングを使って、少ないものは透明にしつつ、どこに点が重なっているかを確認するのも1つかもしれませんし、箱ひげ図で表現するもの良いかもしれません。

g1 <- ggplot(diamonds, aes(x = carat, y = price, col = cut)) +
  geom_point(alpha = 0.1)

g2 <- ggplot(diamonds, aes(x = cut_number(carat, 5), y = price, colour = cut)) +
  geom_boxplot()

grid.arrange(g1, g2, ncol = 2)

練習問題14 : 2次元プロットは、1次元プロットでは見えない異常値を明らかにします。たとえば、下記のプロットでは、xとyの異常な組み合わせがあることがわかります。

アルファブレンディングなので、外れ値でも少ないものは透明に近くなりますが、x:0,y:0の組み合わせは一定数あることがわかります。

ggplot(data = diamonds) +
  geom_point(mapping = aes(x = x, y = y), alpha = 0.2)

5.6 パターンとモデル

練習問題はありません

5.7 ggplot2の呼び出し

練習問題はありません

5.8 さらにまなぶために

練習問題はありません