library(tidyverse)
library(dplyr)
library(nycflights13)
#探索性數(shù)據(jù)分析(exploratory data analysis，EDA)
#? 變量:一種可測量的數(shù)量、質量或屬性。
#? 值:變量在測量時的狀態(tài)。變量值在每次測量之間可以發(fā)生改變。
#? 觀測:或稱個案，指在相同條件下進行的一組測量(通常，一個觀測中的所有測量是在同一時間對同一對象進行的)。
#一個觀測會包含多個值，每個值關聯(lián)到不同的變量。有時我們會將觀測稱為數(shù)據(jù)點。
#? 表格數(shù)據(jù):一組值的集合，其中每個值都關聯(lián)一個變量和一個觀測。
#如果每個值都有自己所屬的“單元”，每個變量都有自己所屬的列，每個觀測都有自己所屬的行，那么表格數(shù)據(jù)就是整潔的。
#變動
#變動是每次測量時數(shù)據(jù)值的變化趨勢。
#所有變量都有自己的變動模式，可以揭示出一些有趣的信息。
#理解這種模式 的最好方法就是對變量值的分布進行可視化表示。
#對分布進行可視化表示
#對變量分布進行可視化的方法取決于變量是分類變量還是連續(xù)變量。
#如果僅在較小的集合內取值，那么這個變量就是分類變量。
#分類變量在 R 中通常保存為因子或字符向量。
#要想檢查分類變量的分布，可以使用條形圖:
ggplot(data = diamonds) +
? geom_bar(mapping = aes(x = cut))
#條形的高度表示每個 x 值中觀測的數(shù)量，你可以使用 dplyr::count() 手動計算出這些值:
diamonds %>%
? count(cut)
#如果可以在無限大的有序集合中任意取值，那么這個變量就是連續(xù)變量。
#數(shù)值型和日期時 間型變量就是連續(xù)變量的兩個例子。
#要想檢查連續(xù)變量的分布，可以使用直方圖:
ggplot(data = diamonds) +
? geom_histogram(mapping = aes(x = carat), binwidth = 0.5)
#你可以通過 dplyr::count() 和 ggplot2::cut_width() 函數(shù)的組合來手動計算結果:
diamonds %>%
? count(cut_width(carat, 0.5))
#直方圖對 x 軸進行等寬分箱，然后使用條形的高度來表示落入每個分箱的觀測的數(shù)量。
#在上圖中，最高的條形表示幾乎有 30 000 個觀測的 carat 值在 0.25 和 0.75 之間，這兩個值分別是條形的左側值和右側值。
#使用 binwidth 參數(shù)來設定直方圖中的間隔的寬度，該參數(shù)是用 x 軸變量的單位來度量的
#在使用直方圖時，你應該試驗一下不同的分箱寬度，因為不同的分箱寬度可以揭示不同的模式。
#如果只考慮重量小于 3 克拉的鉆石，并選擇一個更小的分箱寬度，那么直方圖如下所示:
smaller <- diamonds %>%
? filter(carat < 3)
ggplot(data = smaller, mapping = aes(x = carat)) +
? geom_histogram(binwidth = 0.1)
diamonds %>%
? filter(carat < 3) %>%
ggplot(mapping = aes(x = carat)) +
? geom_histogram(binwidth = 0.1)
#如果想要在同一張圖上疊加多個直方圖，那么我們建議你使用 geom_freqploy() 函數(shù)來代替 geom_histogram() 函數(shù)。
#geom_freqploy() 可以執(zhí)行和 geom_histogram() 同樣的計算過程，
#但前者不使用條形來顯示計數(shù)，而是使用折線。
#疊加的折線遠比疊加的條形更容易理解
ggplot(data = smaller, mapping = aes(x = carat, color = cut)) +
? geom_freqpoly(binwidth = 0.1)
diamonds %>%
? filter(carat < 3) %>%
? ggplot(mapping = aes(x = carat, color = cut)) +
? geom_freqpoly(binwidth = 0.1)
#典型值
#條形圖和直方圖都用比較高的條形表示變量中的常見值，
#而用比較矮的條形表示變量中不常見的值。
#沒有條形的位置表示數(shù)據(jù)中沒有這樣的值。
#? ? 哪些值是最常見的?為什么?
#? ? 哪些值是非常罕見的?為什么?這和你的預期相符嗎?
#? ? 你能發(fā)現(xiàn)任何異乎尋常的模式嗎?如何解釋?
#? ? 為什么重量為整數(shù)克拉和常見分數(shù)克拉的鉆石更多?
#? ? 為什么位于每個峰值稍偏右的鉆石比稍偏左的鉆石更多?
#? ? 為什么沒有重量超過 3 克拉的鉆石?
ggplot(data = smaller, mapping = aes(x = carat)) +
? geom_histogram(binwidth = 0.01)
#一般來說，相似值聚集形成的簇表示數(shù)據(jù)中存在子組。為了理解子組，我們提出以下問題。
#? ? 每個簇中的觀測是如何相似的?
#? ? 不同簇之間的觀測是如何不相似的?
#? ? 如何解釋或描述各個簇?
#? ? 為什么有些簇的外觀可能具有誤導作用?
#異常值
#異常值是與眾不同的觀測或者是模式之外的數(shù)據(jù)點。
#有時異常值是由于數(shù)據(jù)錄入錯誤而產(chǎn)生的;
#有時異常值則能開辟出一塊重要的新科學領域。
#如果數(shù)據(jù)量比較大，有時很難在直 方圖上發(fā)現(xiàn)異常值。
#例如，查看鉆石數(shù)據(jù)集中 y 軸變量的分布，唯一能表示存在異常值的證據(jù)是，y 軸的取值范圍出奇得寬:
ggplot(diamonds) +
? geom_histogram(mapping = aes(x = y), binwidth = 0.5)
#為了更容易發(fā)現(xiàn)異常值，我們可以使用 coord_cartesian() 函數(shù)將 y 軸靠近 0 的部分放大
ggplot(diamonds) +
? geom_histogram(mapping = aes(x = y), binwidth = 0.5) +
? coord_cartesian(ylim = c(0, 50))
#coord_cartesian() 函數(shù)中有一個用于放大 x 軸的 xlim() 參數(shù)
#我們就可以看出有 3 個異常值，分別位于 0、30 左右和 60 左右。
#我們使用 dplyr 將它們找出來:
unusual <- diamonds %>%
filter(y < 3 | y > 20) %>%
? arrange(y)
unusual
#鉆石的寬度不可能是 0 毫米，因此這些值肯定是錯誤的。
#32 毫米和 59 毫米同樣是令人難以置信的， 這樣的鉆石長度超過 1 英寸(1 英寸 =2.54 厘米)，簡直就是無價之寶!
#使用帶有異常值和不帶異常值的數(shù)據(jù)分別進行分析，是一種良好的做法。
#如果兩次分析的結果差別不大，而你又無法說明為什么會有異常值，那么完全可以用缺失值替代異常值，然后繼續(xù)進行分析。
#但如果兩次分析的結果有顯著差別，那么你就不能在沒有正當理由的情況下丟棄它們。
#你需要弄清出現(xiàn)異常值的原因(如數(shù)據(jù)輸入錯誤)，并在文章中說明丟棄它們的理由。
#缺失值
#如果在數(shù)據(jù)集中發(fā)現(xiàn)異常值，但只想繼續(xù)進行其余的分析工作，那么有 2 種選擇。
#將帶有可疑值的行全部丟棄:
diamonds2 <- diamonds %>%
? filter(between(y, 3, 20))
#我們不建議使用這種方式，因為一個無效測量不代表所有測量都是無效的。
#此外，如果數(shù)據(jù)質量不高，若對每個變量都采取這種做法，那么你最后可能會發(fā)現(xiàn)數(shù)據(jù)已經(jīng)所剩無幾!
#相反，我們建議使用缺失值來代替異常值。
#最簡單的做法就是使用 mutate() 函數(shù)創(chuàng)建 一個新變量來代替原來的變量。
#你可以使用 ifelse() 函數(shù)將異常值替換為 NA:
diamonds2 <- diamonds %>%
? mutate(y = ifelse(y < 3 | y > 20, NA, y))
#ifelse() 函數(shù)有 3 個參數(shù)。第一個參數(shù) test 應該是一個邏輯向量，
#如果 test 為 TRUE，函數(shù)結果就是第二個參數(shù) yes 的值;
#如果 test 為 FALSE，函數(shù)結果就是第三個參數(shù) no 的值。
#和 R 一樣，ggplot2 也遵循不能無視缺失值的原則。
#因為無法明確地繪制出缺失值，
#所以 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)
#有時你會想弄清楚造成有缺失值的觀測和沒有缺失值的觀測間的區(qū)別的原因。
#例如，在 nycflights13::flights 中，dep_time 變量中的缺失值表示航班取消了。
#因此，你應該比較 一下已取消航班和未取消航班的計劃出發(fā)時間。
#可以使用 is.na() 函數(shù)創(chuàng)建一個新變量來 完成這個操作:
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
? )
#相關變動
#如果變動描述的是一個變量內部的行為，那么相關變動描述的就是多個變量之間的行為。
#相關變動是兩個或多個變量以相關的方式共同變化所表現(xiàn)出的趨勢。
#查看相關變動的最好方式是將兩個或多個變量間的關系以可視化的方式表現(xiàn)出來。
#如何進行這種可視化表示同樣取決于相關變量的類型。
#分類變量與連續(xù)變量
#我們經(jīng)常需要探索連續(xù)變量的分布，這種分布按照一個分類變量的值可以分為幾個組，
#就像前面的頻率多邊形圖一樣。
#geom_freqpoly() 的默認外觀不太適合這種比較，因為高度是由計數(shù)給出的。
#這就意味著，如果一組觀測的數(shù)量明顯少于其他組的話，就很難看出形狀上的差別。
#舉個例子，我們探索一下鉆石價格是如何隨著質量而變化的:
ggplot(data = diamonds, mapping = aes(x = price)) +
? geom_freqpoly(mapping = aes(color = cut), binwidth = 500)
#很難看出分布上的差別，因為總體看來各組數(shù)量的差別太大了
ggplot(diamonds) +
? geom_bar(mapping = aes(x = cut))
#為了讓比較變得更容易，需要改變 y 軸的顯示內容，不再顯示計數(shù)，而是顯示密度。
#密度是對計數(shù)的標準化，這樣每個頻率多邊形下邊的面積都是 1:
ggplot(
? data = diamonds,
? mapping = aes(x = price, y = ..density..) )+
? geom_freqpoly(mapping = aes(color = cut), binwidth = 500)
#這張圖的部分內容非常令人驚訝，其顯示出一般鉆石(質量最差)的平均價格是最高的!
#但這可能是因為頻率多邊形圖很難解釋，所以這張圖還有很多可以改進的地方。
#按分類變量的分組顯示連續(xù)變量分布的另一種方式是使用箱線圖。
#箱線圖是對變量值分布 的一種簡單可視化表示，這種圖在統(tǒng)計學家中非常流行。
#每張箱線圖都包括以下內容。
#? 一個長方形箱子，下面的邊表示分布的第 25 個百分位數(shù)，上面的邊表示分布的第 75 個 百分位數(shù)，上下兩邊的距離稱為四分位距。
#? 箱子的中部有一條橫線，表示分布的中位數(shù)， 也就是分布的第 50 個百分位數(shù)。
#? 這三條線可以表示分布的分散情況，還可以幫助我們明確數(shù)據(jù)是關于中位數(shù)對稱的，還是偏向某一側。
#? 圓點表示落在箱子上下兩邊 1.5 倍四分位距外的觀測，這些離群點就是異常值，因此需 要單獨繪出。
#? 從箱子上下兩邊延伸出的直線(或稱為須)可以到達分布中最遠的非離群點處。
#使用 geom_boxplot() 函數(shù)查看按切割質量分類的價格分布:
ggplot(data = diamonds, mapping = aes(x = cut, y = price)) +
? geom_boxplot()
#雖然看不出太多關于分布的信息，但箱線圖更加緊湊，因此可以更容易地比較多個類別
#與前面的圖形一樣，我們可以從箱線圖中發(fā)現(xiàn)違反直覺的 現(xiàn)象:質量更好的鉆石的平均價格更低!
#cut 是一個有序因子:“一般”不如“較好”、“較好”不如“很好”，以此類推。
#因為很多分類變量并沒有這種內在的順序，所以有時需要對其重新排序來繪制信息更豐富的圖形。
#重新排序的其中一種方法是使用 reorder() 函數(shù)。
#我們看一下 mpg 數(shù)據(jù)集中的 class 變量。
#你可能很想知道公路里程因汽車類別的不 同會有怎樣的變化:
ggplot(data = mpg, mapping = aes(x = class, y = hwy)) +
? geom_boxplot()
#為了更容易發(fā)現(xiàn)趨勢，可以基于 hwy 值的中位數(shù)對 class 進行重新排序:
ggplot(data = mpg) +
? geom_boxplot(
? ? mapping = aes(
? ? ? x = reorder(class, hwy, FUN = median),
? ? ? y = hwy
? ? ) )
#如果變量名很長，那么將圖形旋轉 90 度效果會更好一些。
#你可以通過 coord_flip() 函數(shù)完成這一操作:
ggplot(data = mpg) +
? geom_boxplot(
? ? mapping = aes(
? ? ? x = reorder(class, hwy, FUN = median),
? ? ? y = hwy
? ? ) )+
? coord_flip()
#兩個分類變量
#要想對兩個分類變量間的相關變動進行可視化表示，
#需要計算出每個變量組合中的觀測數(shù)量。
#完成這個任務的其中一種方法是使用內置的 geom_count() 函數(shù)
ggplot(data = diamonds) +
? geom_count(mapping = aes(x = cut, y = color))
#圖中每個圓點的大小表示每個變量組合中的觀測數(shù)量。
#相關變動就表示為特定 x 軸變量值 與特定 y 軸變量值之間的強相關關系。
#計算變量組合中的觀測數(shù)量的另一種方法是使用 dplyr:
diamonds %>%
? count(color, cut)
#接著使用 geom_tile() 函數(shù)和填充圖形屬性進行可視化表示:
diamonds %>%
? count(color, cut) %>%
? ggplot(mapping = aes(x = color, y = cut)) +
? geom_tile(mapping = aes(fill = n))
#如果分類變量是無序的，那么可以使用 seriation 包對行和列同時進行重新排序，以便更清楚地表示出有趣的模式。
#對于更大的圖形，你可以使用 d3heatmap 或 heatmaply 包，這兩個包都可以生成有交互功能的圖形。
# 一般會將大數(shù)量的分類變量放在y軸或者名字比較長的放y軸
#兩個連續(xù)變量
#對于兩個連續(xù)變量間的相關變動的可視化表示，我們已經(jīng)介紹了一種非常好的方法:
#使用 geom_point() 畫出散點圖。你可以將相關變動看作點的模式。
#例如，你可以看到鉆石的克 拉數(shù)和價值之間存在一種指數(shù)關系:
ggplot(data = diamonds) +
? geom_point(mapping = aes(x = carat, y = price))
#隨著數(shù)據(jù)集規(guī)模的不斷增加，散點圖的用處越來越小，
#因為數(shù)據(jù)點開始出現(xiàn)過繪制，并堆積在一片黑色區(qū)域中(如上面的散點圖所示)。
#使用 alpha 圖形屬性添加透明度
ggplot(data = diamonds) +
? geom_point(
? ? mapping = aes(x = carat, y = price),
? ? alpha = 1 / 100
? )
#但是很難對特別大的數(shù)據(jù)集使用透明度。
#另一種解決方法是使用分箱。
#之前使用了 geom_histogram() 和 geom_freqpoly() 函數(shù)在一個維度上進行分箱
#geom_bin2d() 和 geom_hex() 函數(shù)在兩個維度上進行分箱
#geom_bin2d() 和 geom_hex() 函數(shù)將坐標平面分為二維分箱，
#并使用一種填充顏色表示落入每個分箱的數(shù)據(jù)點。
#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ù)變量的作用就相當于分類變量。
#例如，你 可以對 carat 進行分箱，然后為每個組生成一個箱線圖:
ggplot(data = smaller, mapping = aes(x = carat, y = price)) +
? geom_boxplot(mapping = aes(group = cut_width(carat, 0.1)))
#以上示例使用了 cut_width(x, width) 函數(shù)將 x 變量分成寬度為 width 的分箱。
#默認情況下，不管其中有多少個觀測，箱線圖看上去都差不多(除了離群點的數(shù)量不同)
#因此很 難分辨出每個箱線圖是對不同數(shù)量的觀測進行摘要統(tǒng)計的。
#如果想要體現(xiàn)這種信息，可以使用參數(shù) varwidth = TRUE 讓箱線圖的寬度與觀測數(shù)量成正比。
ggplot(data = smaller, mapping = aes(x = carat, y = price)) +
? geom_boxplot(mapping = aes(group = cut_width(carat, 0.1)),varwidth = TRUE)
#另一種方法是近似地顯示每個分箱中的數(shù)據(jù)點的數(shù)量，此時可以使用 cut_number() 函數(shù)
ggplot(data = smaller, mapping = aes(x = carat, y = price)) +
? geom_boxplot(mapping = aes(group = cut_number(carat, 20)))
#二維圖形可以顯示一維圖形中看不到的離群點。
#例如，以下圖形中的有些點具有異常的 x 值和 y 值組合，這使得這些點成為了離群點，
#即使這些點的 x 值和 y 值在單獨檢驗時似乎是正常的。
ggplot(data = diamonds) +
? geom_point(mapping = aes(x = x, y = y)) +
? coord_cartesian(xlim = c(4, 11), ylim = c(4, 11))
#模式和模型
#數(shù)據(jù)中的模式提供了關系線索。
#如果兩個變量之間存在系統(tǒng)性的關系，那么這種關系就會在數(shù)據(jù)中表現(xiàn)為一種模式。
#模式是數(shù)據(jù)科學中最有效的工具之一，因為其可以揭示相關變動。
#如果說變動會生成不確定性，那么相關變動就是減少不確定性。
#如果兩個變量是共同變化的，就可以使用一個變 量的值來更好地預測另一個變量的值。
#如果相關變動可以歸因于一種因果關系(一種特殊 情況)，那么就可以使用一個變量的值來控制另一個變量的值。
#模型是用于從數(shù)據(jù)中抽取模式的一種工具。
#例如，我們思考一下鉆石數(shù)據(jù)。
#切割質量與價格之間的關系是很難理解的，因為切割質量和克拉數(shù)以及克拉數(shù)和價格之間是緊密相關的。
#我們可以使用模型去除價格和克拉數(shù)之間的強關系，這樣就可以繼續(xù)研究剩余的微妙關系。
#以下代碼擬合了一個模型，可以根據(jù) carat 預測 price，并計算出殘差(預測值和 實際值之間的差別)。
#一旦去除克拉數(shù)對價格的影響，殘差就能反映出鉆石的價格:
library(modelr)
#擬合了一個模型，可以根據(jù) carat 預測 price
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))
#去除克拉數(shù)和價格之間的強關系后，就可以看到預料中的切割質量與價格的關系，
ggplot(data = diamonds2) +
? geom_boxplot(mapping = aes(x = cut, y = resid))
#對于同樣大小的鉆石，切割質量更好的鉆石更昂貴:
#ggplot2調用
#介紹 ggplot2 代碼的一種更精簡表示
ggplot(data = faithful, mapping = aes(x = eruptions)) +
? geom_freqpoly(binwidth = 0.25)
#一個函數(shù)的前一個或前兩個參數(shù)是非常重要的，你應該將它們牢記于心。
#ggplot() 函數(shù)的前兩個參數(shù)是 data 和 mapping，aes() 函數(shù)的前兩個參數(shù)是 x 和 y。
#在本書剩余的部分中，我們不再寫出這些參數(shù)名，
#這樣既可以節(jié)省輸入時間，也可以讓代碼樣板更精簡，以便更容易找出兩張圖之間的不同之處。
#以更精簡的方式重寫前面的繪圖語句，結果如下所示:
ggplot(faithful, aes(eruptions)) +
? geom_freqpoly(binwidth = 0.25)
#有時我們會將數(shù)據(jù)轉換管道操作的最終結果轉換為一張圖。
#注意，此時轉換是用 + 號實現(xiàn)的，不是 %>%。
#我也不希望如此，但遺憾的是，ggplot2 是在管道操作方式發(fā)明前開發(fā)出來的。
diamonds %>%
? count(cut, clarity) %>%
? ggplot(aes(clarity, cut, fill = n)) +
? geom_tile()
