R for Data Scienceの例題を解く- Chapter 7 探索的データ分析

神Hadley R for Data Science の例題たちとその解答を書き残します。
今回はChapter 7 Explatory Data Analysisです。

他の章の例題と解答たちはこちら

この章ではライブラリ群tidyverseを使います。

library(tidyverse)

Chapter 7 Explatory Data Analysis

7.3 Variation

7.3.4 Exercises

1. Explore the distribution of each of the x, y, and z variables in diamonds. What do you learn? Think about a diamond and how you might decide which dimension is the length, width, and depth.

ggplot(diamonds) + geom_histogram(aes(x = x), bins = 30)
ggplot(diamonds) + geom_histogram(aes(x = y), bins = 30)
ggplot(diamonds) + geom_histogram(aes(x = z), bins = 30)

xのヒストグラム

yのヒストグラム

zのヒストグラム

y, zでは少数ではあるが極端に値が大きな測定値があることが分かります。
縦にトンガッているダイアモンドは考えにくいので、外れ値の無いxが高さではないか。
横幅と奥行きとでは横幅の方が伸ばしやすいような気がするので、外れ値がより大きなyが横幅、残ったzが奥行きだと推測しました。


2. Explore the distribution of price. Do you discover anything unusual or surprising? (Hint: Carefully think about the binwidth and make sure you try a wide range of values.)

 ggplot(diamonds) +
   geom_histogram(aes(price), bins = 300) +
   scale_x_continuous(minor_breaks = seq(0, 20000, by = 1000))
priceのヒストグラム

価格の分布はゼロから$500~$1000周辺にあるピークに向けて急上昇し、そこからなだらかに減少していきます。
ピークは最もポピュラーなダイアモンドの価格を示していると考えられます。

diamonds %>%
  count(cut_number(price, 300)) %>%
  filter(n == max(n))
## A tibble: 1 x 2
#  `cut_number(price, 300)`     n
#                    <fctr> <int>
#1                (770,776]   261

$770あたりが普通のダイアモンドと言えそう。

この分布の形は自然であり納得できますが、いくつかの例外も見いだせます。

  1. $1500あたりに分布のギャップがある。
  2. $1000~$3000の間にピークがぴょんぴょんと立っている。
  3. $4000を過ぎたあたりで二度目の山が立っている

まず1点目のギャップは最も目立つ例外です。
具体的にどの値なのかを調べます。

 ggplot(diamonds) +
   geom_histogram(aes(price), bins = 500) +
   coord_cartesian(x = c(1000, 2000))
ダイアモンド価格の分布図をギャップ周辺で拡大

$1500ドル周辺にはっきりとしたギャップがあることが観察できます。
$1500ドルで売ってはならないという協定があるのか、$1500ドルのダイヤモンドは縁起が悪いというジンクスがあるのか、そうではなく実は$1500周辺の値だけがデータ収集の過程で欠落してしまったのではないかなどということが考えられます。


3. How many diamonds are 0.99 carat? How many are 1 carat? What do you think is the cause of the difference?

diamonds %>% summarise(sum(carat == 0.99), sum(carat == 1))
## A tibble: 1 x 2
#  `sum(carat == 0.99)` `sum(carat == 1)`
#                 <int>             <int>
#1                   23              1558

1カラットという数字は大台であるため、それを上回るかどうかはダイアモンドの価値を大きく動かすと考えられます。

diamonds %>%
  group_by(carat) %>%
  summarise(mean_price = mean(price)) %>%
  ggplot(aes(carat, mean_price)) +
    geom_point() +
    coord_cartesian(x = c(0.9,1.1), y = c(2500,7500))
ダイアモンドのカラットと平均価格の関係図

このように価格に及ぼす影響が大きいので、カットで削る量を調整してでも1カラットを下回らないようにしているのではないか。


4. Compare and contrast coord_cartesian() vs xlim() or ylim() when zooming in on a histogram. What happens if you leave binwidth unset? What happens if you try and zoom so only half a bar shows?

coord_cartesianはズームアップするだけだが、xlim, ylimはグラフを描画する対象を制限するという違いがあります。
だからcoord_cartesianbinwidthに影響を与えないが、xlimは指定した範囲内でヒストグラムの幅を設定しようとします。
またylimでy座標の範囲を指定すると、その範囲から漏れたデータは描画されなくなるので、ヒストグラムの棒が立ちません。
一方coord_cartesianではズームアップしているだけなのでヒストグラムの棒は立ちます。


7.4 Missing values

7.4.1 Exercises

What happens to missing values in a histogram? What happens to missing values in a bar chart? Why is there a difference?

nycflights13::flights %>% ggplot(aes(arr_delay)) + geom_histogram(binwidth = 1)
nycflights13::flights %>% ggplot(aes(arr_delay)) + geom_histogram(binwidth = 1, na.rm = TRUE)
nycflights13::flights %>% ggplot(aes(arr_delay)) + geom_bar()
nycflights13::flights %>% ggplot(aes(arr_delay)) + geom_bar(na.rm = TRUE)

私の目にはこれらの4つのチャートは全く同じに見えますが、何が違うのでしょうか?
NAは単にプロットされないという点で同じように思います。


2. What does na.rm = TRUE do in mean() and sum()?

NAは除外してNA以外で平均、及び合計を計算します。

7.5 Covariation

7.5.1 A categorical and continuous variable

7.5.1.1 Exercises

1. Use what you’ve learned to improve the visualisation of the departure times of cancelled vs. non-cancelled flights.

nycflights13::flights %>%
   mutate(
        dep_min = dep_time %/% 100 * 60 + dep_time %% 100,
        cancelled = is.na(arr_time)
        ) %>%
   ggplot(aes(dep_min, ..density..)) +
     geom_freqpoly(aes(color = cancelled), na.rm = TRUE, binwidth = 30) +
     scale_x_continuous(minor_breaks=seq(0, 1500, by = 100))
フライトが欠航したかどうかごとの出発時刻の分布図。横軸の単位はAM0時からの経過分

上記のグラフより一日の前半は欠航便の確率密度の方が少ないが、15時よりも前あたりで逆転し夕方以降は欠航便の確率密度の方が逆転します。
すなわち午後では便が欠航する確率の方が高いです。


2. What variable in the diamonds dataset is most important for predicting the price of a diamond? How is that variable correlated with cut? Why does the combination of those two relationships lead to lower quality diamonds being more expensive?

diamonds %>% ggplot(aes(carat, price)) + geom_point()
ダイアモンドの大きさと価格との関係図。正の相関が見て取れる。

大きさ(carat)が価格と相関が大きいです。
ではcaratcutとの相関はどうでしょうか。

diamonds %>% ggplot(aes(cut, carat)) + geom_boxplot() + coord_flip()
ダイアモンドのカットごとの価格の箱ひげ図

cutfairの場合が顕著にcaratが大きいことが分かります。
このためfaircaratに引っ張られて価格が大きくなっていると推測できます。


3. Install the ggstance package, and create a horizontal boxplot. How does this compare to using coord_flip()?

install.packages("ggstance")
diamonds %>% ggplot(aes(carat, cut)) + ggstance::geom_boxploth()
diamonds %>% ggplot(aes(cut, carat)) + geom_boxplot()

coord_flipと違うのはmappingxyとが逆転していること。


4. One problem with boxplots is that they were developed in an era of much smaller datasets and tend to display a prohibitively large number of “outlying values”. One approach to remedy this problem is the letter value plot. Install the lvplot package, and try using geom_lv() to display the distribution of price vs cut. What do you learn? How do you interpret the plots?

install.packages("lvplot")
diamonds %>% ggplot(aes(cut, carat)) + geom_lv() + coord_flip()
ダイアモンドの大きさと価格の分布をLV箱グラフで図示

LV箱グラフを見るとcutVery Goodの価格分布が最も足が短いことが分かります。
通常の箱ひげ図ではGoodVery Goodとのどちらが分布の足が短いのか見出すことは難しいです。
このように分布の足の長さを比較するためには箱ひげ図よりもLV箱グラフの方が優れています。


5. Compare and contrast geom_violin() with a facetted geom_histogram(), or a coloured geom_freqpoly(). What are the pros and cons of each method?

diamonds %>%
   mutate(cut = reorder(cut, as.numeric(cut) * -1)) %>%
   ggplot(aes(cut, carat)) +
     geom_violin() +
     coord_flip()
diamonds %>%
   ggplot(aes(carat, ..density..)) +
     geom_histogram(bins = 20) +
     facet_wrap(~cut, ncol = 1, strip.position = "left") +
     scale_y_continuous(minor_breaks = NULL) +
     theme(strip.text.y = element_text(angle = 180))
diamonds %>%
   ggplot(aes(carat, ..density..)) +
     geom_freqpoly(aes(color = cut)) +
     theme(legend.position = "left")

cutとcaratとのviolinプロット

caratのヒストグラムをcutで分割

caratのfreqpolyをcut別に表示

method pro con
バイオリンプロット 裾の長さを比較しやすい ある値における分布の大小を比較しにくい。
分割ヒストグラム グラフの軸を個別に設定できる。 分布の比較をしにくい
freqpoly 分布の大小の比較をしやすい 個々の分布の形がわかりにくい

6. If you have a small dataset, it’s sometimes useful to use geom_jitter() to see the relationship between a continuous and categorical variable. The ggbeeswarm package provides a number of methods similar to geom_jitter(). List them and briefly describe what each one does.

install.packages("ggbeeswarm")
library(ggbeeswarm)
ggplot(mpg, aes(class, hwy)) + geom_quasirandom()
ggplot(mpg, aes(class, hwy)) + geom_jitter()

geom_quasirandommpgclass, hwyをプロット

geom_jittermpgclass, hwyをプロット

  • geom_quasirandomは実現値が(ほぼ)左右対称に分布していますが、geom_jitterでは実現値に対称性がなく分布はランダム
  • geom_quasirandomはブレが左右だけであるが、geom_jitterでは上下左右に散布しています。ただしgeom_jitterのオプションとしてheight = 0とすれば左右のみの散布になります。

7.5.2 Two categorical variables

7.5.2.1 Exercises

1. How could you rescale the count dataset above to more clearly show the distribution of cut within colour, or colour within cut?

# 元々のcount plot
diamonds %>%
  count(color, cut) %>%  
  ggplot(mapping = aes(x = color, y = cut)) +
    geom_tile(mapping = aes(fill = n))

# 分布にエッジをかけたcount plot
diamonds %>%
  count(color, cut) %>%
  ggplot(mapping = aes(x = color, y = cut)) +
    geom_tile(mapping = aes(fill = n^2))

ダイヤモンドの色とカットのカウントプロット。明るい方がデータ数が多い。

ダイヤモンドの色とカットのカウントプロット。明るい方がデータ数が多い。色のコントラストを極端にした。

2. Use geom_tile() together with dplyr to explore how average flight delays vary by destination and month of year. What makes the plot difficult to read? How could you improve it?

nycflights13::flights %>%
   group_by(dest, month) %>%
   summarise(arr_delay = mean(arr_delay, na.rm = TRUE)) %>%
   ggplot(aes(dest, factor(month))) +
     geom_tile(aes(fill = arr_delay)) +
     theme(axis.text.x = element_text(angle = 90))
flightsの行き先と月別で平均遅れ時間をプロット。遅れが大きいほど色が明るい

データに隙間があってすべての色が塗られていないのがデータを見えにくくしているように思えます。
よってデータを寄せてみます。

nycflights13::flights %>%
   mutate(dest = reorder(dest, arr_delay, sum, na.rm = TRUE)) %>%
   group_by(dest, month) %>%
   summarise(arr_delay = mean(arr_delay, na.rm = TRUE)) %>%
   ggplot(aes(dest, factor(month))) +
     geom_tile(aes(fill = arr_delay)) +
     theme(axis.text.x = element_text(angle = 90))

flightsの行き先と月別で平均遅れ時間をプロット。遅れが大きいほど色が明るい。行き先の表示順を遅れが多い順に並び替えた.。

これをみれば1月、7月、8月が周辺の月よりも色が明るい、すなわち遅れが発生している、ということがわかります。
横軸もCAE(Columbia Metolopolitan Airport)の成績の悪さも目立ちます。


3. Why is it slightly better to use aes(x = color, y = cut) rather than aes(x = cut, y = color) in the example above?

colorは7種類、cutは5種類であり、colorの方が多い。
一般的にコンピュータのディスプレイは横長なので種類数が多いほうを横軸に置くほうが良い。


7.5.3 Two continuous variables

7.5.3.1 Exercises

1. Instead of summarising the conditional distribution with a boxplot, you could use a frequency polygon. What do you need to consider when using cut_width() vs cut_number()? How does that impact a visualisation of the 2d distribution of carat and price?

frequency polygonはヒストグラムのように分布の広がりを表します。
ヒストグラムを使うとき、比較したいグループに含まれるサンプル数が同じでなければ分布を比較できません。
そのためにはcut_numberを使います。

smaller <- diamonds %>% filter(carat < 3) 
smaller %>%  ggplot(aes(price)) + geom_freqpoly(aes(color = cut_number(carat, 10)), bins = 30)
ダイヤモンドの価格の分布。カラット別。ただしそれぞれのカラット集合には同じ数のサンプルが入っている

これを見ればcaratが大きくなるほど価格の分布が広がっていく、という傾向が簡単に読み取れます。

また比較対象を等間隔で刻みたいという場合は、縦軸をカウントではなく密度にすればスケールは合うようにできます。

smaller %>% ggplot(aes(price, ..density..)) + geom_freqpoly(aes(color = cut_width(carat, 0.5)))
ダイヤモンドの価格の分布。カラット別。ただしそれぞれのカラット集合は等幅で区切っている

これでも比較することはできますが、最大のグループはサンプル数が3個しかないため意味のある比較が難しい。


2. Visualise the distribution of carat, partitioned by price.

diamonds %>%
   ggplot(aes(carat)) +
     geom_freqpoly(aes(color = cut_number(price, 5)), bins = 20)
ダイヤモンドのカラットの分布。価格別

それぞれの価格帯ごとにピークとなる大きさがあることが分かる。
ただし最高級クラスだけは1, 1.5, 2カラットの3つの山を持っている。


3. How does the price distribution of very large diamonds compare to small diamonds. Is it as you expect, or does it surprise you?

diamonds %>%
  ggplot(aes(cut_number(carat, 8), price)) +
    geom_violin() + coord_flip()
ダイヤモンドの価格の分布。カラット別

一番上のグループはカラットの範囲が1.34-5.01と広いので分布が広がるのは理解できます。
しかし上から二番目のグループはカラットの範囲は1.04-1.34とあまり広くはないにも関わらず一番上のグループとあまり変わらない分布です。
まるで1カラットよりも大きくなると価格と大きさの相関が無くなってしまうようです。

またダイヤモンドがいくら大きくなっても19000ドルは絶対に越えられない壁として立ちはだかっています。
1カラットでも、3カラットでも同じです。

diamonds %>%
   ggplot(aes(carat, price)) +
     geom_point(aes(color = cut_number(carat, 6))) +
     geom_hline(yintercept = 19000, linetype = 2, size = 1)
ダイヤモンドのカラットと価格の散布図。大きくなっても19000ドルは越えられない

4. Combine two of the techniques you’ve learned to visualise the combined distribution of cut, carat, and price.

diamonds %>%
  ggplot(aes(cut_number(carat, 8), price)) +
    geom_boxplot(aes(fill = cut), size = 0.3, outlier.size = 0.5) +
    theme(axis.text.x = element_text(angle = 90))
ダイヤモンドの価格の箱ひげ図。カラットとカット別

5. Two dimensional plots reveal outliers that are not visible in one dimensional plots. For example, some points in the plot below have an unusual combination of x and y values, which makes the points outliers even though their x and y values appear normal when examined separately.

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

分布図にするということは一変数で見るということ、つまりデータを一次元に射影するということです。
それに対して散布図はデータを二次元に射影している。
このため外れ値を見つけることができます。
三次元プロットを使えばもっと緻密に外れ値を見つけられるはずです。


関連記事

参考サイト

R for Data Science Chapter 7 Exploratory Data Analysis

コメントを残す