明確概念：探索性數(shù)據(jù)分析（exploratory data analysis， EDA），一般過程為：
(1) 對數(shù)據(jù)提出問題。
(2) 對數(shù)據(jù)進行可視化、轉(zhuǎn)換和建模，進而找出問題的答案。
(3) 使用上一個步驟的結(jié)果來精煉問題，并提出新問題。

5.3.1　對分布進行可視化表示

確定變量是分類變量還是連續(xù)變量，要想檢查分類變量的分布，可以使用條形圖：

library(tidyverse)
head(diamonds)
ggplot(diamonds,aes(x=cut))+ geom_bar() 

條形的高度表示每個 x 值中觀測的數(shù)量，可以使用 dplyr::count() 手動計算出這些值：

diamonds%>% count(cut)

要想檢查連續(xù)變量的分布，可以使用直方圖：

ggplot(data = diamonds) +
geom_histogram(mapping = aes(x = carat), binwidth = 0.5)

可以通過 dplyr::count() 和 ggplot2::cut_width() 函數(shù)的組合來手動計算結(jié)果. binwidth 參數(shù)來設定直方圖中的間隔的寬度，該參數(shù)是用 x 軸變量的單位來度量的。

diamonds %>%
count(cut_width(carat, 0.5))

smaller <- diamonds %>% filter(carat < 3)

ggplot(diamonds,aes(carat)) + geom_histogram(bindwidth=0.1)

在同一張圖上疊加多個直方圖, 用geom_freqploy()代替geom_histogram()，用折線表示。

ggplot(smaller,aes(carat,color=cut)) + geom_freqpoly(binwidth=0.1)

相似值聚集形成的簇表示數(shù)據(jù)中存在子組。

5.3.3　異常值

ggplot(diamonds) + geom_histogram(aes(y),bindwidth=0.5)
# 可以看出y軸有一大片空白，因為有異常值出現(xiàn)
# 放大y軸看一下
ggplot(diamonds) + 
  geom_histogram(aes(y),binwidth = 0.5) +
  coord_cartesian(ylim = c(0,50))

# 用dplyr中的filter將異常值找出來
unusual <- diamonds %>% filter(y<3 | y>20) %>% arrange(y)
unusual

coord_cartesian() 函數(shù)中有一個用于放大 x 軸的 xlim() 參數(shù)。 ggplot2 中也有功能稍有區(qū)
別的 xlim() 和 ylim() 函數(shù)：它們會忽略溢出坐標軸范圍的那些數(shù)據(jù)。
如果帶有異常值和不帶異常值的數(shù)據(jù)分別進行分析，結(jié)果差別較大的話要找出異常值的原因，如果差別不大，可以用NA代替。

5.3.4　練習
(1)研究 x、 y 和 z 變量在 diamonds 數(shù)據(jù)集中的分布。你能發(fā)現(xiàn)什么？思考一下，對于一條
鉆石數(shù)據(jù)，如何確定表示長、寬和高的變量？

ggplot(diamonds) + geom_density(aes(z))
head(diamonds)
diamonds %>%
  mutate(id = row_number()) %>%
  select(x, y, z, id) %>%
  gather(variable, value, -id)  %>%
  ggplot(aes(x = value)) +
  geom_density() +
  geom_rug() +
  facet_grid(variable ~ .)

(2)研究 price 的分布，你能發(fā)現(xiàn)不尋?；蛄钊梭@奇的事情嗎？（提示：仔細考慮一下
binwidth 參數(shù)，并確定試驗了足夠多的取值。）

head(diamonds)
diamonds  %>% ggplot() + geom_histogram(aes(price),binwidth = 10)
 ggplot(filter(diamonds,price<2500)) + geom_histogram(aes(price),binwidth = 10)
 
 ggplot(filter(diamonds), aes(x = price)) +
  geom_histogram(binwidth = 100, center = 0)
 
 diamonds %>%
  mutate(ending = price %% 10) %>%
  ggplot(aes(x = ending)) +
  geom_histogram(binwidth = 1, center = 0) +
  geom_bar()
 
 diamonds %>%
  mutate(ending = price %% 100) %>%
  ggplot(aes(x = ending)) +
  geom_histogram(binwidth = 1) +
  geom_bar()
 
 diamonds %>% mutate(ending=price%%1000) %>% filter(ending>500, ending<1000) %>%
   ggplot(aes(x=ending)) +
   geom_histogram(binwidth = 1) + geom_bar()

(3) 0.99 克拉的鉆石有多少？ 1 克拉的鉆石有多少？造成這種區(qū)別的原因是什么？

head(diamonds)
diamonds %>% filter(carat >= 0.99, carat <= 1) %>% count(carat)
# 發(fā)現(xiàn)0.99的比較少，是因為0.01的差別，1克拉賣的比較貴嗎？可以看下其他類型的數(shù)量
diamonds %>% filter(carat >= 0.9, carat <= 1.1) %>% count(carat) %>% print(n=30)

(4)比較并對比 coord_cartesina() 和 xlim()/ylim() 在放大直方圖時的功能。如果不設置
binwidth 參數(shù)，會發(fā)生什么情況？如果將直方圖放大到只顯示一半的條形，那么又會發(fā)
生什么情況？

ggplot(diamonds) + geom_histogram(aes(price)) + coord_cartesian(xlim = c(100,5000),ylim = c(0,3000))
# 下面這種方法將不在xlim和ylim范圍的值都去掉了，而coord_cartesian方法是計算后取值，并沒有將范圍外的點去掉。
ggplot(diamonds) + geom_histogram(aes(price)) + xlim(100,5000)+ylim(0,3000)

5.4　缺失值

數(shù)據(jù)中有異常值，可以將異常值去掉：

dim(diamonds)
diamonds2 <- diamonds %>% filter(between(y,3,20))
dim(diamonds2)

一般不建議去掉，建議使用缺失值來代替異常值。

diamonds2 <- diamonds %>% mutate(y=ifelse(y<3 | y>20, NA, y))

ifelse函數(shù)參數(shù)1放入邏輯判斷，如果為T，結(jié)果就是第二個參數(shù)的值，如果為F，就是第三個參數(shù)的值。
ggplot2會忽略缺失值：

ggplot(diamonds2, aes(x,y)) + geom_point()
# na.rm=TRUE 可以去掉缺失值
# is.na將含NA的行進行區(qū)分
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
)

練習
(1) 直方圖如何處理缺失值？條形圖如何處理缺失值？為什么會有這種區(qū)別？

diamonds2 <- diamonds %>%
  mutate(y = ifelse(y < 3 | y > 20, NA, y))

ggplot(diamonds2, aes(x = y)) +
  geom_histogram()
# 直方圖對x進行計數(shù)，自動去掉缺失值
# 柱狀圖將NA也統(tǒng)計為一個變量進行計數(shù)
diamonds %>%
  mutate(cut = if_else(runif(n()) < 0.1, NA_character_, as.character(cut))) %>%
  ggplot() +
  geom_bar(mapping = aes(x = cut))

(2) na.rm = TRUE 在 mean() 和 sum() 函數(shù)中的作用是什么？
移除缺失值再進行統(tǒng)計

5.5　相關變動

5.5.1　分類變量與連續(xù)變量

# geom_freqpoly對于變異較小的變量不容易看出分別
# 這里頻率多邊形圖顯示的是數(shù)量
ggplot(diamonds,aes(price)) + geom_freqpoly(aes(color=cut),binwidth=500)

# 可以將y軸轉(zhuǎn)變?yōu)槊芏?ggplot(diamonds,aes(price,..density..)) + geom_freqpoly(aes(color=cut),binwidth=500)

按分類變量的分組顯示連續(xù)變量分布的另一種方式是使用箱線圖

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

ggplot(mpg,aes(class,hwy)) + geom_boxplot()
# 為了更容易發(fā)現(xiàn)趨勢，可以基于 hwy 值的中位數(shù)對 class 進行重新排序
ggplot(mpg,aes(x=reorder(class,hwy,FUN=median),hwy)) + geom_boxplot()

# coord_flip將圖形旋轉(zhuǎn) 90 度
ggplot(mpg,aes(x=reorder(class,hwy,FUN=median),hwy)) + geom_boxplot() + coord_flip()

練習
(1) 前面對比了已取消航班和未取消航班的出發(fā)時間，使用學習到的知識對這個對比的可視
化結(jié)果進行改善。

# 比較取消和未取消的航班出發(fā)時間，簡化為分類變量與連續(xù)變量的關系，用boxplot
head(nycflights13::flights)
nycflights13::flights %>% mutate( canceled = 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() +
    geom_boxplot(mapping = aes(y = sched_dep_time, x = canceled))

(2) 在鉆石數(shù)據(jù)集中，哪個變量對于預測鉆石的價格最重要？這個變量與切割質(zhì)量的關系是
怎樣的？為什么這兩個變量的關系組合會導致質(zhì)量更差的鉆石價格更高呢？

# 很明顯carat與price的價格最相關了，這兩個都是連續(xù)變量，可以用點圖表示
ggplot(diamonds,aes(x=carat,y=price)) + geom_point()

# 也可以將x分割為一個個單元進行統(tǒng)計
ggplot(diamonds,aes(carat,price)) + geom_boxplot(aes(group=cut_width(carat,0.1)))

# 查看顏色與價格的關系
ggplot(diamonds,aes(color,price)) + geom_boxplot()

ggplot(data = diamonds) +
  geom_boxplot(mapping = aes(x = clarity, y = price))

(3) 安裝 ggstance 包，并創(chuàng)建一個橫向箱線圖。這種方法與使用 coord_flip() 函數(shù)有何區(qū)別？

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

library("ggstance")
# 可以看出x和y軸進行了調(diào)換
ggplot(data = mpg) +
  geom_boxploth(mapping = aes(y = reorder(class, hwy, FUN = median), x = hwy))

(4) 箱線圖存在的問題是，在小數(shù)據(jù)集時代開發(fā)而成，對于現(xiàn)在的大數(shù)據(jù)集會顯示出數(shù)量極
其龐大的異常值。解決這個問題的一種方法是使用字母價值圖。安裝 lvplot 包，并嘗試
使用 geom_lv() 函數(shù)來顯示價格基于切割質(zhì)量的分布。你能發(fā)現(xiàn)什么問題？如何解釋這
種圖形？

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

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

(5) 比較并對比 geom_violin()、分面的 geom_histogram() 和著色的 geom_freqploy()。每種方法的優(yōu)缺點是什么？

ggplot(data = diamonds, mapping = aes(x = price, y = ..density..)) +
  geom_freqpoly(mapping = aes(color = cut), binwidth = 500)

ggplot(diamonds,aes(price)) + geom_histogram() + facet_wrap(~ cut, ncol = 1, scales = "free_y")

ggplot(data = diamonds, mapping = aes(x = cut, y = price)) +
  geom_violin() +
  coord_flip()

(6) 對于小數(shù)據(jù)集，如果要觀察連續(xù)變量和分類變量間的關系，有時使用 geom_jitter() 函數(shù)是特別有用的。 ggbeeswarm 包提供了和 geom_jitter()相似的一些方法。列出這些方法
并簡單描述每種方法的作用。

library("ggbeeswarm")
ggplot(data = mpg) +
  geom_quasirandom(mapping = aes(x = reorder(class, hwy, FUN = median),
                                 y = hwy))

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

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

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

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

5.5.2　兩個分類變量

兩個分類變量的關系肯定要先計數(shù)，可以用geom_count()函數(shù)
d3heatmap 或 heatmaply 包可以生成交互式圖

ggplot(diamonds) + geom_count(aes(cut,color))
str(diamonds)
# dplyr也可以用來計數(shù)
diamonds %>% count(color,cut)
#geom_tile函數(shù)填充熱圖
diamonds %>% count(color,cut) %>% ggplot(aes(color,cut)) +
  geom_tile(aes(fill=n))

練習
(1) 如何調(diào)整 count 數(shù)據(jù)，使其能更清楚地表示出切割質(zhì)量在顏色間的分布，或者顏色在切
割質(zhì)量間的分布？

# 另外一種表示方法就是求比例
library(viridis)
diamonds %>% count(color,cut) %>% group_by(color) %>% mutate(prop=n/sum(n)) %>% ggplot(aes(color,cut)) + geom_tile(aes(fill=prop)) + scale_fill_viridis(limits=c(0,1))

diamonds %>%
  count(color, cut) %>%
  group_by(cut) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(mapping = aes(x = color, y = cut)) +
  geom_tile(mapping = aes(fill = prop)) +
  scale_fill_viridis(limits = c(0, 1))

(2) 使用 geom_tile() 函數(shù)結(jié)合 dplyr來探索平均航班延誤數(shù)量是如何隨著目的地和月份的
變化而變化的。為什么這張圖難以閱讀？如何改進？

nycflights13::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() +
  labs(x = "Month", y = "Destination", fill = "Departure Delay")
# 可以看出信息很雜亂，可以對目的地排序，移除缺失值，改變顏色等進行改進
nycflights13::flights    %>% group_by(month,dest) %>% summarise(dep_delay=mean(dep_delay,na.rm=T)) %>% 
  group_by(dest) %>% filter(n() == 12) %>% ungroup() %>% mutate(dest=reorder(dest,dep_delay)) %>% 
  ggplot(aes(x=factor(month),y=dest,fill=dep_delay)) +
  geom_tile() + scale_fill_viridis() + labs(x = "Month", y = "Destination", fill = "Departure Delay")

(3) 為什么在以上示例中使用 aes(x = color, y = cut) 要比 aes(x = cut, y = color) 更好？

# 一般會將大數(shù)量的分類變量放在y軸或者名字比較長的放y軸
diamonds %>%
  count(color, cut) %>%  
  ggplot(mapping = aes(y = color, x = cut)) +
    geom_tile(mapping = aes(fill = n))

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

5.5.3　兩個連續(xù)變量

連續(xù)變量之間的關系一般用散點圖來表示。geom_point()

ggplot(data = diamonds) +
geom_point(mapping = aes(x = carat, y = price))
# 
ggplot(data = diamonds) +
geom_point(
mapping = aes(x = carat, y = price),
alpha = 1 / 100
)

對于大數(shù)據(jù)集，為了避免重合，可以用geom_bin2d() 和 geom_hex()函數(shù)將坐標平面分為二維分箱，并使用一種填充顏色表示落入
每個分箱的數(shù)據(jù)點。

ggplot(diamonds,aes(carat,price)) + geom_bin2d()

另一種方式是對一個連續(xù)變量進行分箱，因此這個連續(xù)變量的作用就相當于分類變量。
cut_width(x, width) 函數(shù)將 x 變量分成寬度為 width 的分箱。參數(shù) varwidth = TRUE 讓箱線圖的寬度與觀測數(shù)量成正比。
cut_number() 函數(shù)近似地顯示每個分箱中的數(shù)據(jù)點的數(shù)量

ggplot(diamonds,aes(carat,price)) +geom_boxplot(aes(group=cut_width(carat,0.1))) 

ggplot(diamonds,aes(carat,price)) +geom_boxplot(aes(group=cut_number(carat,20))) 

練習
(1) 除了使用箱線圖對條件分布進行摘要統(tǒng)計，你還可以使用頻率多邊形圖。使用 cut_
width() 函數(shù)或 cut_number() 函數(shù)時需要考慮什么問題？這對 carat 和 price 的二維分
布的可視化表示有什么影響？

ggplot(diamonds,aes(color=cut_width(carat,1,boundary = 0),x=price)) + 
  geom_freqpoly() +ylab('Carat')

(2) 按照 price 分類對 carat 的分布進行可視化表示。

ggplot(diamonds,aes(cut_number(price,10),y=carat)) + geom_boxplot() +coord_flip() 

(3) 比較特別大的鉆石和比較小的鉆石的價格分布。結(jié)果符合預期嗎？還是出乎意料？

(4) 組合使用你學習到的兩種技術(shù)，對 cut、 carat 和 price 的組合分布進行可視化表示。

(5) 二維圖形可以顯示一維圖形中看不到的離群點。例如，以下圖形中的有些點具有異常的
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))

5.6　模式和模型

數(shù)據(jù)中的模式提供了關系線索,用于探索兩個變量的相關性。
模型是用于從數(shù)據(jù)中抽取模式的一種工具。
殘差（預測值和實際值之間的差別）

library(modelr)
mod <- lm(log(price) ~ log(carat),data=diamonds)

diamonds2 <- diamonds %>% add_residuals(mod) %>% mutate(resid=exp(resid))

ggplot(diamonds2) + geom_point(aes(carat,resid))

閱讀推薦：
生信技能樹公益視頻合輯：學習順序是linux，r，軟件安裝，geo，小技巧，ngs組學！
B站鏈接：https://m.bilibili.com/space/338686099
YouTube鏈接：https://m.youtube.com/channel/UC67sImqK7V8tSWHMG8azIVA/playlists
生信工程師入門最佳指南：https://mp.weixin.qq.com/s/vaX4ttaLIa19MefD86WfUA