該部分學習內(nèi)容來自《R for Data Science》。
這次我們學習可視化與處理數(shù)據(jù)來系統(tǒng)地探索數(shù)據(jù)——統(tǒng)計學家稱之為探索性數(shù)據(jù)分析（exploratory data analysis），簡稱為EDA。
EDA是一個迭代的圓圈：
- 生成關于你所擁有數(shù)據(jù)的問題
- 通過對數(shù)據(jù)可視化、轉(zhuǎn)換和建模尋找答案
- 使用你學到的重定義你的問題或者生成新的問題
相比于嚴格的規(guī)則與流程，EDA更像一種思考狀態(tài)。在EDA的初始階段，你可以隨意探索跳入你腦海的任意一個想法。
數(shù)據(jù)清理僅僅是EDA的一個應用，你詢問是否你的數(shù)據(jù)如你所期。想要進行數(shù)據(jù)清理，你將需要使用EDA所有的工具：可視化、轉(zhuǎn)換與建模。
準備
這部分我們將使用dplyr與ggplot2交互地對數(shù)據(jù)提問并回答，請確保在學習前安裝好tidyverse包，可以順利進行以下操作：
library(tidyverse)
## -- Attaching packages ------------------------------------------------------------------ tidyverse 1.2.1 --
## <U+221A> ggplot2 2.2.1 <U+221A> purrr 0.2.4
## <U+221A> tibble 1.4.2 <U+221A> dplyr 0.7.4
## <U+221A> tidyr 0.8.0 <U+221A> stringr 1.3.0
## <U+221A> readr 1.1.1 <U+221A> forcats 0.3.0
## -- Conflicts --------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
EDA的目標是理解數(shù)據(jù)，最容易實現(xiàn)的方式是用問題作為用具引導我們進行探索。想要提出有質(zhì)量的問題首先必須提出大量問題。
提問的方式?jīng)]有準則可言，但兩類基本問題的提出會對我們理解數(shù)據(jù)大有裨益：
- 在我的變量中有什么類型的變異？
- 在我的變量之間有什么共同的變異？
在正式進行分析之前，允許我定義一些術語：
- 一個變量是一個你可以測量的數(shù)量、質(zhì)量或者屬性
- 一個值是當你測量一個變量時它的狀態(tài)
- 一個觀測或者一個案例是相似條件下一組測量值的集合。一次觀測通常包含多個值，每個值與不同的變量相關聯(lián)。我將有時稱一個觀測為一個數(shù)據(jù)點
- 表格數(shù)據(jù)是值的集合，每一個值都與一個變量和一個觀測相關聯(lián)
可視化分布
你如何可視化一個變量的分布取決于該變量是連續(xù)還是分類的。如果一個變量僅能取少量的幾個值，我們就說它是分類變量。在R中，分類變量常保存為因子或者是字符向量。想要檢查分類變量的分布，使用直方圖：
ggplot(data = diamonds) +
geom_bar(mapping = aes(x = cut))
直方圖的高度顯示了每一個x值有多少個觀測。你可以使用dplyr::count()手動進行計算：
diamonds %>%
count(cut)
## # A tibble: 5 x 2
## cut n
## <ord> <int>
## 1 Fair 1610
## 2 Good 4906
## 3 Very Good 12082
## 4 Premium 13791
## 5 Ideal 21551
或者
diamonds %>%
group_by(cut) %>%
summarise(n = n())
## # A tibble: 5 x 2
## cut n
## <ord> <int>
## 1 Fair 1610
## 2 Good 4906
## 3 Very Good 12082
## 4 Premium 13791
## 5 Ideal 21551
想要檢查連續(xù)變量的分布，使用直方圖：
ggplot(data = diamonds) +
geom_histogram(mapping = aes(x=carat), binwidth = 0.5)
你也可以結合dplyr::count()與ggplot2::cut_width()手動計算這個：
diamonds %>%
count(cut_width(carat, 0.5))
## # A tibble: 11 x 2
## `cut_width(carat, 0.5)` n
## <fct> <int>
## 1 [-0.25,0.25] 785
## 2 (0.25,0.75] 29498
## 3 (0.75,1.25] 15977
## 4 (1.25,1.75] 5313
## 5 (1.75,2.25] 2002
## 6 (2.25,2.75] 322
## 7 (2.75,3.25] 32
## 8 (3.25,3.75] 5
## 9 (3.75,4.25] 4
## 10 (4.25,4.75] 1
## 11 (4.75,5.25] 1
直方圖將x軸劃分為相等的寬度（由bin控制），然后使用每個條形的高度來代表落入該區(qū)域觀測的數(shù)目。條形的寬度非常重要，不同的設定可能揭示出數(shù)據(jù)的內(nèi)在分布模式。
smaller <- diamonds %>%
filter(carat < 3)
ggplot(data = smaller, mapping = aes(x = carat)) + geom_histogram(binwidth = 0.1)
如果你想要在同一幅圖上畫多個直方圖，我推薦你使用geom_freqpoly()函數(shù)而不是geom_histogram()，它們的內(nèi)部計算一致，但前者使用線形。
ggplot(data = smaller, mapping = aes(x = carat, color = cut)) +
geom_freqpoly(binwidth = 0.1)
典型值
在條形圖和直方圖中，高的條形顯示了變量通常出現(xiàn)的值，而短的條形說明比較少見，而沒有條形出現(xiàn)的地方說明變量幾乎不可能會出現(xiàn)那個值。想要把這些信息轉(zhuǎn)為有用的問題，尋找不同尋常之處：
- 哪個值最常出現(xiàn)？為什么?
- 哪個值很少出現(xiàn)？為什么？這符合你的預期嗎?
- 你可以看到任何異常的模式嗎?什么可以解釋它？
作為一個示例，下面的直方圖表明一些有趣的問題：
- 為什么在整個克拉和普通的克拉中有更多的鉆石？
- 為什么每個峰頂右側的鉆石略多于每個峰頂左側的鉆石？
- 為什么沒有鉆石大于3克拉？
ggplot(data = smaller, mapping = aes(x = carat)) +
geom_histogram(binwidth = 0.01)
通常，有相似值的集群表明你數(shù)據(jù)中存在亞組（可以分組）。想要理解亞組，提問：
- 每個集群內(nèi)的觀察結果如何相似？
- 如何在不同的集群中觀察到彼此不同的結果？
- 你如何解釋或描述集群？
- 為什么集群的外觀會產(chǎn)生誤導？
下面的直方圖顯示了黃石國家公園老忠實噴泉272次噴發(fā)的長度（以分鐘為單位）。噴發(fā)時間似乎分為兩組：噴發(fā)時間短（約2分鐘）和噴發(fā)時間長（4-5分鐘），但間隔很少。
ggplot(data = faithful, mapping = aes(x = eruptions)) +
geom_histogram(binwidth = 0.25)
上面的許多問題都會提示你探索變量之間的關系，例如，查看一個變量的值是否可以解釋另一個變量的行為。 我們很快會做到。
不尋常的值
異常值是不尋常的觀察結果; 數(shù)據(jù)點似乎不符合模式。有時異常值是數(shù)據(jù)輸入錯誤;其他時間異常值表明重要的新科學發(fā)現(xiàn)。 當你有很多數(shù)據(jù)時，在直方圖中有時很難看到異常值。例如，從鉆石數(shù)據(jù)集中分配y變量。異常值的唯一證據(jù)是x軸上異常寬的限制。
ggplot(diamonds) +
geom_histogram(mapping = aes(x = y), binwidth = 0.5)
在常見的箱子有很多觀察結果，罕見的箱子都很短，以至于看不到它們（盡管也許如果你專注于0，你會發(fā)現(xiàn)一些東西）。 為了便于查看異常值，我們需要使用coord_cartesian()縮放y軸的小值：
ggplot(diamonds) +
geom_histogram(mapping = aes(x = y), binwidth = 0.5) +
coord_cartesian(ylim = c(0, 50))
（coord_cartesian()同樣也有一個xlim參數(shù)當你在需要縮小x軸時可以使用。ggplot2也有xlim()與ylim()函數(shù)，但工作不太一樣：它們會扔掉所有超過該限制的所有數(shù)據(jù)。）
這讓我們可以看到存在3個異常值：一個是0，一個30左右，一個60左右。我們使用dplyr把它們?nèi)〕觯?/p>
unusual <- diamonds %>%
filter(y < 3 | y > 20) %>%
select(price, x, y, z) %>%
arrange(y)
unusual
## # A tibble: 9 x 4
## price x y z
## <int> <dbl> <dbl> <dbl>
## 1 5139 0. 0. 0.
## 2 6381 0. 0. 0.
## 3 12800 0. 0. 0.
## 4 15686 0. 0. 0.
## 5 18034 0. 0. 0.
## 6 2130 0. 0. 0.
## 7 2130 0. 0. 0.
## 8 2075 5.15 31.8 5.12
## 9 12210 8.09 58.9 8.06
y變量測量這些鉆石的三個維度之一，單位為mm。 我們知道鉆石不能有0mm的寬度，所以這些值肯定是不正確的。 我們也可能懷疑32毫米和59毫米的尺寸是不合理的：那些鉆石長1英寸，但價格沒有數(shù)十萬美元！
無論是否存在異常值，重復進行分析都是很好的做法。如果他們對結果的影響最小，并且您無法弄清楚他們?yōu)槭裁磿霈F(xiàn)這種情況，那么將其替換為缺失值并繼續(xù)前進是合理的。但是，如果它們對您的結果有重大影響，則沒有理由不放棄它們。你需要弄清楚是什么導致了他們（例如數(shù)據(jù)輸入錯誤），并且你在寫作中應透露刪除了它們。
練習
- 探索鉆石中每個x，y和z變量的分布。你學到什么？想想鉆石，以及哪個尺寸是長度，寬度和深度。
- 探索價格的分布。你發(fā)現(xiàn)什么不尋?；蛄钊梭@訝？（提示：仔細考慮binwidth，并確保您嘗試各種值。）
- 0.99克拉多少顆鉆石？ 1克拉多少錢？ 你認為什么是差異的原因？
- 在放大直方圖時比較和對比
coord_cartesian()與xlim()或ylim()。 如果您保留binwidth未設置會發(fā)生什么？ 如果您嘗試縮放以便只顯示一半，會發(fā)生什么情況？
# 1
# x 分布
ggplot(data = diamonds) + geom_histogram(mapping = aes(x = x), binwidth = 0.5)
ggplot(data = diamonds) + geom_histogram(mapping = aes(x = y), binwidth = 0.5)
ggplot(data = diamonds) + geom_histogram(mapping = aes(x = z), binwidth = 0.5)
# 發(fā)現(xiàn)它們的分布模式其實非常接近
# 鉆石長寬深有什么特點我不是很清楚，只能按常理x,y,z分別對應長寬深
# 2 探索價格分布
ggplot(data = diamonds) + geom_histogram(mapping = aes(x = price))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# 基本趨勢是價格越高的數(shù)量越少，除了x軸5000前面有點反常
# 3
diamonds %>% filter(carat == 0.99) %>% count
## # A tibble: 1 x 1
## n
## <int>
## 1 23
diamonds %>% filter(carat == 0.99) %>% select(price, carat) %>%
ggplot(aes(x=price)) + geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
diamonds %>% filter(carat == 1) %>% select(price, carat) %>%
ggplot(aes(x=price)) + geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# 4 我們以書中的例子做比較
ggplot(diamonds) +
geom_histogram(mapping = aes(x = y), binwidth = 0.5) +
coord_cartesian(ylim = c(0, 50))
ggplot(diamonds) +
geom_histogram(mapping = aes(x = y)) + ylim(0, 50)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 3 rows containing missing values (geom_bar).
缺失值
如果你發(fā)現(xiàn)數(shù)據(jù)集中存在缺失值，想要簡便地進行后續(xù)的分析，你有兩個選項：
- 丟掉出現(xiàn)異常值的整行
diamonds2 <- diamonds %>%
filter(between(y, 3, 20))
我并不推薦這樣操作，因為一個測量值出現(xiàn)問題并不代表所有數(shù)據(jù)都有問題。另外，如果你的數(shù)據(jù)質(zhì)量比較低，可能會刪除絕大多數(shù)的數(shù)據(jù)，那你還分析什么呢？
- 相比于上述方法，我推薦用缺失值替換異常值。最簡單的方式是使用
mutate()去創(chuàng)建一個修改的版本，你可以使用ifelse()函數(shù)將NA替換異常值。
diamonds2 <- diamonds %>%
mutate(y = ifelse(y < 3| y >20, NA, y))
ifelse()函數(shù)有三個參數(shù)，第一個參數(shù)需要是邏輯向量，當test某個值為TRUE時，對應結果就是第二個參數(shù)，反之是第三個參數(shù)。
像R，ggplot2應用了缺失值從不該在沉默中缺失的哲學。當你繪圖時很難發(fā)現(xiàn)缺失值，所以ggplot2會提醒你那么被移除了。
ggplot(data = diamonds2, mapping = aes(x=x, y=y))+
geom_point()
## Warning: Removed 9 rows containing missing values (geom_point).
想要抑制警告，設置na.rm=TRUE：
ggplot(data = diamonds2, mapping = aes(x=x, y=y)) +
geom_point(na.rm = TRUE)
有時候你會想什么讓觀測缺失了數(shù)值，它與有記錄的值又有什么不同。例如，在nycflights13::flights數(shù)據(jù)集中，dep_time變量的缺失值（參加使用dplyr數(shù)據(jù)處理一文）顯示了航班被取消了。所以你想比較安排起飛取消和沒取消的次數(shù)：
nycflights13::flights %>%
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(mapping = aes(sched_dep_time)) +
geom_freqpoly(mapping = aes(color = cancelled), binwidth = 1/4)
共變化
如果變化描述變量內(nèi)的行為，則共變化描述變量之間的行為。 共變是兩個或更多變量的值以相關方式一起變化的趨勢。 發(fā)現(xiàn)共變化的最好方法是可視化兩個或更多變量之間的關系。 你怎么做，應該再次取決于所涉及的變量的類型。
一個分類變量和連續(xù)變量
如前一個頻率多邊形一樣，想要探索由分類變量細分的連續(xù)變量的分布是很常見的。geom_freqpoly()的默認外觀對于那種比較沒有用處，因為高度由count給出。這意味著如果其中一個組別比其他組別小得多，很難看出形狀上的差異。 例如，讓我們來探討鉆石的價格如何隨其質(zhì)量而變化：
ggplot(data = diamonds, mapping = aes(x = price)) +
geom_freqpoly(mapping = aes(color=cut), binwidth = 500)
It’s hard to see the difference in distribution because the overall counts differ so much:
ggplot(diamonds) +
geom_bar(mapping = aes(x = cut))
為了使比較更容易，我們需要交換y軸上顯示的內(nèi)容。 我們不顯示計數(shù)，而是顯示密度，這是標準化的計數(shù)，以便每個頻率多邊形下的面積為1。
ggplot(data = diamonds, mapping = aes(x = price,
y = ..density..)) + geom_freqpoly(mapping = aes(color = cut), binwidth = 500)
這個結果有一些令人驚訝的地方 -看起來，一般的鉆石（最低的質(zhì)量）的平均價格最高！ 但也許這是因為頻率多邊形有點難以解釋 - 在這個結果中有很多事情要做。
顯示分類變量細分的連續(xù)變量分布的另一種方法是boxplot。 boxplot是一種在統(tǒng)計學家中很受歡迎的價值分布的視覺縮寫。 每個boxplot包括：
- 一個從分布的第25個百分點延伸到第75個百分點的盒子，這個距離被稱為四分位間距（IQR）。在框的中間是顯示分布的中值（即第50百分位）的線。這三條線讓你了解分布的擴散以及分布是否對稱于中值或偏向一側。
- 顯示觀察結果的如果點從盒子的任一邊緣落入IQR的1.5倍以上。 這些偏離點是不尋常的，因此單獨繪制。
- 從盒子兩端延伸出來的線（或胡須）分配到最遠的非異常點。
讓我們使用geom_boxplot()來看一下價格的分布：
ggplot(data = diamonds, mapping = aes(x=cut, y=price)) +
geom_boxplot()
我們看到關于分布的信息要少得多，但箱形圖更加緊湊，因此我們可以更輕松地比較它們（并且更適合于一個圖）。 它支持反直覺的發(fā)現(xiàn)，即更好質(zhì)量的鉆石平均價格更便宜！
許多分類變量沒有這樣的內(nèi)在順序，因此您可能需要對它們重新排序以提供更豐富的信息。一種方法是使用reorder()函數(shù)。
例如，使用mpg數(shù)據(jù)集中的class變量。您可能有興趣了解不同類別汽車的公路里程如何變化：
ggplot(data = mpg, mapping = aes(x=class, y=hwy)) +
geom_boxplot()
想要使趨勢更清楚，我們可以使用hwy的中位數(shù)進行排序：
ggplot(data = mpg) +
geom_boxplot(mapping = aes(x = reorder(class, hwy, FUN = median), y=hwy))
如果你由一個長的變量名，你反轉(zhuǎn)90度顯示的效果更佳：
ggplot(data = mpg) +
geom_boxplot(mapping = aes(x = reorder(class, hwy, FUN=median), y=hwy)) + coord_flip()
兩個分類變量
想要可視化兩個分類變量的共變化，你需要對每一種組合進行計數(shù)。一種方法是依賴內(nèi)置的geom_count()：
ggplot(data = diamonds) +
geom_count(mapping = aes(x=cut, y=color))
圖中每個圓圈的大小顯示每個值組合處發(fā)生的觀察次數(shù)。 共變將表現(xiàn)為具體x值與具體y值之間的強相關性。
另外一種方法是使用dplyr計算計數(shù)：
diamonds %>%
count(color, cut)
## # A tibble: 35 x 3
## color cut n
## <ord> <ord> <int>
## 1 D Fair 163
## 2 D Good 662
## 3 D Very Good 1513
## 4 D Premium 1603
## 5 D Ideal 2834
## 6 E Fair 224
## 7 E Good 933
## 8 E Very Good 2400
## 9 E Premium 2337
## 10 E Ideal 3903
## # ... with 25 more rows
使用geom_tile()和fill映射可視化：
diamonds %>%
count(color, cut) %>%
ggplot(mapping = aes(x=color, y=cut)) +
geom_tile(mapping = aes(fill = n))
如果分類變量無序，則可能需要使用seriation包同時對行和列進行重新排序，以便更清楚地顯示有趣的模式。 對于較大的地塊，你可能需要嘗試使用d3heatmap或heatmaply創(chuàng)建交互式地的圖形。
兩個連續(xù)變量
你已經(jīng)見過一種可以非常好可視化兩個連續(xù)變量的方式——用geom_point()畫一個散點圖。你可以通過這些點來看變量共變化的模式。例如你可以看到鉆石的大小與其價格有指數(shù)關系。
ggplot(data = diamonds) +
geom_point(mapping = aes(x=carat, y=price))
當數(shù)據(jù)集越來越大時，散點圖會變得越來越?jīng)]有，因為點會覆蓋糾纏到一塊。一種優(yōu)化的方式是使用alpha映射來添加透明度。
ggplot(data = diamonds) +
geom_point(mapping = aes(x=carat, y=price), alpha = 1 / 100)
但數(shù)據(jù)集如果太大效果也不好。另一個辦法是使用geom_bin2d()和geom_hex()對二維數(shù)據(jù)分組。
這兩個函數(shù)可以從二維角度將某一塊區(qū)域坐標的點劃分為一塊，使用顏色等映射進行填充可視化。geom_bin2d()創(chuàng)建矩形，而geom_hex()創(chuàng)建六邊形。在使用geom_hex()之前需要先安裝hexbin包。
ggplot(data = smaller) +
geom_bin2d(mapping = aes(x = carat, y = price))
#install.packages("hexbin")
ggplot(data = smaller) +
geom_hex(mapping = aes(x = carat, y = price))
還有一種處理兩個連續(xù)變量的方式是將一個連續(xù)變量（進行分組）作為分類變量看待，然后使用我們前面提到的可視化一個分類變量與一個連續(xù)變量的方式。例如：
ggplot(data = smaller, mapping = aes(x=carat, y=price)) +
geom_boxplot(mapping = aes(group = cut_width(carat, 0.1)))
上面使用的cut_width(x, width)將x根據(jù)寬度width劃分了許多組別。默認不管有多少個觀測值箱線圖看起來都一樣（除了離群點），所以很難說明每個箱子匯總了多少數(shù)據(jù)點。我們可以將箱子的（寬度）大小和匯總觀測值的數(shù)目成比例，設置varwidth=TRUE即可。
另外一種方式是盡量讓每個箱子觀察值都差不多，這可以通過cut_number()實現(xiàn)：
ggplot(data = smaller, mapping = aes(x=carat, y=price)) +
geom_boxplot(mapping = aes(group = cut_width(carat, 0.1)), varwidth = TRUE)
ggplot(data = smaller, mapping = aes(x = carat, y=price)) + geom_boxplot(mapping = aes(group = cut_number(carat, 20)))
模式與模型
數(shù)據(jù)的模式提供了（變量）關系線索。如果兩個變量存在系統(tǒng)性的關系，它必定會在數(shù)據(jù)中顯示一定的模式。如果你想要找到一個模式，請這一問自己：
- 這個模式的出現(xiàn)是偶然嗎？
- 你怎么描述這個模式顯示的關系？
- 這個模式顯示的關系有多強？
- 有其他變量可能會影響這個關系嗎？
- 如果你只單獨看一小組的數(shù)據(jù)，這個關系會改變嗎？
老忠實噴泉等待時間與噴發(fā)顯示了這樣一個模式：等待時間越長噴發(fā)也就越長。散點圖可以顯示我們上面注意到的兩個簇團：
ggplot(data = faithful) +
geom_point(mapping = aes(x = eruptions, y= waiting))
模型是從數(shù)據(jù)中抽取出模式的工具。下面代碼擬合了一個鉆石大小carat預測價格price并計算殘差的模型。
library(modelr)
mod <- lm(log(price) ~ log(carat), data=diamonds)
diamonds2 <- diamonds %>%
add_residuals(mod) %>%
mutate(resid = exp(resid))
ggplot(data = diamonds2) +
geom_point(mapping = aes(x = carat, y = resid))
一旦你移除這種強相關關系，你會發(fā)現(xiàn)鉆石質(zhì)量越好價格越貴：
ggplot(data = diamonds2) +
geom_boxplot(mapping = aes(x = cut, y = resid))
Related
