社會科學統計方法

資料視覺化

1 課程目標

我們將介紹R的繪圖指令，包括ggplot2的套件中的指令，把資料視覺化。雖然之前我們已經看過很多例子，但是仍然有必要詳細介紹包括直方圖、折線圖、散佈圖等等圖形，結合dplyr套件中的指令，以培養視覺化資料的能力，並且在圖形上面顯示統計結果。例如我們可以視覺化資料的分佈並且找出每一個類別的平均數：

file<-here("data","studentsfull.txt")
library(ggplot2)
students<-read.table(file, sep='', header=TRUE)
p <- ggplot(students, aes(x=Department, y=Score)) +
  geom_point() +
  theme_classic() 

p + ggplot2::stat_summary(fun.y='mean', colour='red',size=3, geom = "point")

---

2 基礎繪圖指令

有人說，「一張照片勝過千言萬語」。越來越來媒體報導、研究報告、論文以圖形表現變數之間的關係，讓讀者可以很快地了解。

R內建了許多基本的繪圖指令，讓我們可以很輕鬆地使用，但是ggplot2提升了繪圖的層次，所以我們必須學習並且熟練ggplot2。

2.1 基本元素

在R的基礎繪圖功能中，繪圖的元素有「點」、「線」、「大小」、「粗細」以及「顏色」等等。

• cex: 控制符號的大小
• col: 控制符號的顏色，也可以是灰階的濃度
• pch: 控制符號的形狀
• lwd: 控制線的粗細
• lty: 控制線的形狀

圖 2.1表示10種不同濃度黑色、10種不同形狀以及2種大小的點狀的散佈圖，加上粗細、顏色、形狀不同的垂直與水平線。

xn<-c(1:10)
yn<-c(1:10)
col<-c(80, 75, 70, 60, 50, 40, 30, 20, 10, 0)
graycol<-paste("gray",col,sep="")
graycol

##  [1] "gray80" "gray75" "gray70" "gray60" "gray50" "gray40" "gray30" "gray20"
##  [9] "gray10" "gray0"

plot(xn, yn,
  cex=c(rep(1, 5), rep(2,5)),
  col=graycol,
  pch=c(1,2,5,6,8, 16, 17, 20, 22, 24),
  xlab="pch", ylab="cex",
  xaxt="n", yaxt="n",
  xlim=c(0,11),ylim=c(0,11), font.lab=2, cex.lab=1.5)
axis(1, at=c(1:10), labels=c(1,2,5,6,8, 16, 17, 20, 22, 24))
axis(2, at=c(1:10), labels=c(rep(1, 5), rep(2,5)))
abline(v=6, lwd=2, col='red'); abline(h=6, lwd=3, lty=2)

Figure 2.1: 圖型的元素

從以上的指令以及圖形不難看出一些常用圖形參數。另外要說明幾個設定參數：

• xlab: X軸的標題
• ylab: Y軸的標題
• main: 圖形標題
• xlim: X軸最高的值
• ylim: Y軸最高的值
• xaxt='n': 控制X軸不顯示任何數字
• yaxt='n': 控制Y軸不顯示任何數字
• axis: X軸或是Y軸顯示的數字（labels ）以及位置（at ）
• font.lab: 控制X, Y軸名稱的大小
• cex.lab: 控制X, Y軸刻度名稱的大小

那麼我們如何應用到資料分析？例如我們有一筆資料anscombe，其中有兩個連續變數x1, y1，我們可以用散佈圖 2.2 來檢視：

with(anscombe, plot(x1, y1, pch=16, cex=2, col='blue',
     xlab='X', ylab='Y', main='x1-y1', xaxt='n'))
axis(side=1, labels=c(1,"" ,"" ,"" ,"" ,6), at=c(4,6,8,10,12,14))

Figure 2.2: 散佈圖一

在第一個步驟中，圖 2.2的X軸沒有刻度，在第二步驟，我們加上頭跟尾兩個刻度。換句話說，透過axis()函數，我們可以設定X或Y軸任意的刻度。

我們也可以先幫觀察值畫一個框，再用points這個函數，有同樣的效果，如圖2.3：

with(anscombe, summary(x1))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     4.0     6.5     9.0     9.0    11.5    14.0

with(anscombe, summary(y1))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    4.26    6.32    7.58    7.50    8.57   10.84

x<-c(1: max(anscombe$x1))
y<-c(1: max(anscombe$x1))
plot(x, y, type='n')
with(anscombe, points(x1, y1, pch=16, cex=2, col='gray25'))

Figure 2.3: 散佈圖二

我們可以用points函式繪圖標示兩個變數對應的數值，並且用which()控制散佈圖。 圖2.4顯示的顏色或者形狀，類似加上第三個變數的效果：

x<-c(1: max(anscombe$x1))
y<-c(1: max(anscombe$x1))
plot(x, y, type='n')

with(anscombe, points(x1[which(x1<=8)], 
        y1[which(x1<=8)], pch=16, cex=2, col='red'))
with(anscombe, points(x1[which(x1>8)], 
      y1[which(x1>8)], pch=10, cex=2, col='#ee0022'))
abline(v=8, lwd=2, lty=3)

Figure 2.4: 散佈圖三

圖 2.4 顯示x1超過8是加上十字的圓點，小於8是紅色的點。

• 我們可以用colors()這個函式顯示657種顏色。
• 另外我們可以用#加上6碼的Hex指定#RRGGBB的顏色。可以查這個網站輸入Hex6位數找到喜歡的顏色。
• 有時候我們想增加多一點資訊，可以用text在觀察點旁邊加上文字，讓讀者更容易瞭解觀察值的順位。pos 控制文字的位置。例如圖 2.5 顯示11個觀察值的落點。

plot(x, y, type='n')
with(anscombe, points(x1[which(x1<=8)], 
        y1[which(x1<=8)], pch=16, cex=2, col='red'))
with(anscombe, points(x1[which(x1>8)], 
        y1[which(x1>8)], pch=16, cex=2, col='darkblue'))
abline(v=8, lwd=2, lty=3)
with(anscombe, text(x1, y1, 
            c(1: nrow(anscombe)), pos=4))



Figure 2.5: 散佈圖四

• 也可以貼上觀察值本身的數字如圖 2.6：

plot(x, y, type='n', xlim=c(1, 16), xaxt='n')
axis(1, labels = c(2:15), at=c(2:15))
with(anscombe, points(x1[which(x1<=8)], y1[which(x1<=8)], pch=16, 
                  cex=2, col='red'))
with(anscombe, points(x1[which(x1>8)], y1[which(x1>8)], pch=22, 
                cex=2, col='darkblue'))
abline(v=8, lwd=2, lty=3)
with(anscombe, text(x1, y1, paste(x1, y1, sep=","), pos=4))



Figure 2.6: 散佈圖五

☛ 請根據ISLR::College這筆資料的Top10perc排序，篩選出新生為高中前10%畢業生的比率最高的前10名學校，觀察錄取率(Accept/Apps)與曾捐款的校友比率(perc.alumni)的關係，然後註記學校的名稱。

---

2.2 直方圖

直方圖可以表現單一類別變數的分佈，讓人一眼就可以看出哪一個類別有最多的次數或者百分比。例如圖 2.7 顯示學生的科系分佈：

TXT<-here("data","studentsfull.txt")
students<-read.table(TXT, header=T, sep='')
stu <- table(students$Department)
barplot(stu, main="Departments", 
xlab="", ylab="frequency")

Figure 2.7: 直方圖一

直方圖的對象應該是具有名稱的摘要，所以要先把數字、字串或者是類別的向量改為表格的形式。直方圖也可以指定直方本身的顏色以及外框的顏色，例如圖 2.8 換了新的顏色。

barplot(100*stu/nrow(students), main="Departments", 
      xlab="", ylab="Percent", border='red', 
      col='#0011EE22', ylim=c(0, 20))

Figure 2.8: 直方圖二

如果想要表現兩個變數的交集，例如經濟系裡面有多少男生或女生的比例，可以先產生一個交叉列表，然後用prop.table()的函數產生條件機率，就可以繪出堆疊的直方圖 2.9。注意我們加上了圖例，才能了解不同色塊代表的意義。

如果圖例所使用的變數是類別變數，圖例函數寫成： legend = levels(students$Gender))，也就是回傳該變數的類別。如果在讀這筆資料時設定保留字串變數、不轉成類別變數。在此就要把性別轉為類別變數，才能帶出兩個類別，也就是legend = levels(as.factor(students$Gender))。

student.table <- table(students$Gender, students$Department)
barplot(100*prop.table(student.table, margin=2), 
         col=c('brown', 'white'),
        legend = levels(students$Gender))

Figure 2.9: 直方圖三

我們可以換一個圖例的製作方式。圖 2.10 先畫直方圖， 再另外加上圖例，函數內需要設定文字及顏色，但是不能設定符號。請注意圖例的顏色與文字必須與圖形一致。

圖例的位置有“bottomright”, “bottom”, “bottomleft”, “left”, “topleft”, “top”, “topright”, “right”, “center”等等，比較常用的是topleft, topright, bottomright, bottomleft四種。

student.table <- table(students$Gender, students$Department)
barplot(100*prop.table(student.table, margin=2), 
         col=c('blue', '#EE330011'))
legend("top", fill=c('blue', '#EE330011'), c("M","F"))

Figure 2.10: 直方圖與圖例

• 我們也可以設定圖的外框：par()，把圖例放在圖的外面：

par(xpd=T, mar=par()$mar+c(0,0,0,6))
TXT<-here::here("data","studentsfull.txt")
students<-read.table(TXT, header=T, sep='')
stu <- table(students$Department)
student.table <- table(students$Gender, students$Department)
barplot(student.table, col=c("red", "blue"),
  legend=c('female','male'))

Figure 2.11: 圖例在圖外面

pt.students<-100*prop.table(student.table, 2)
barplot(pt.students, col=c("salmon1", "royalblue1"),
        ylim=c(0,100),xlim=c(0.3,8))
legend(8.5, 90, c('female','male'), 
       col=c("salmon1", "royalblue1"),
       pch=c(20,20), bty='n')

Figure 2.12: 圖例在圖外面

如果想要更改類別的名稱，可以直接在barplot()之中設定參數，例如在直方圖 2.13，我們修改了類別的名稱：

barplot(100*prop.table(student.table, margin=2),
      names.arg=c("Aer","Che.", "Eco.", "Eng.", "Jou.", "Mec.","Phy.")) 

Figure 2.13: 直方圖四

☛請嘗試讀取PP0797B2.sav這筆資料，然後畫長條圖表示partyid這個變數中，政黨各類別的相對次數。

---

2.3 盒型圖 (Box Plot)

盒型圖又稱為箱型圖，於1977年由普林斯頓大學統計系教授約翰·圖基（John Tukey）發明。可以表現單一連續變數的中位數、25分位數、75分位數、極端值等等，也可以比較不同類別之下，連續變數的分佈情形。

例如我們觀察mpg這個變數的盒型圖如圖 2.14：

boxplot(mtcars$mpg, ylim=c(0,40), yaxt='n')

Figure 2.14: 箱型圖一

中間的粗線代表中位數。盒型圖上方的線稱為inner fence或者是whisker，代表數列的極大值與75分位數加上四分位距的1.5倍兩個其中的極小值，下方的線代表數列的極小值與25分位數減去四分位距的1.5倍兩個其中的極大值。四分位距則是75分位數減去25分位數。超出inner fences被稱為極端值（outlier）。

為了確認以上的定義，計算這個變數的25與75分位數：

quantile(mtcars$mpg, c(.25, .5, .75), type=6)

##   25%   50%   75% 
## 15.27 19.20 22.80

mtcars裡面的mpg的IQR約等於7.6，所以75分位數加上四分位距的1.5倍與25分位數減去1.5倍四分位距可用R計算如下：

qu<-quantile(mtcars$mpg, c(.25, .5, .75), type=7)
qu

##   25%   50%   75% 
## 15.43 19.20 22.80

upper<- qu[3]+1.5*(qu[3]-qu[1])
lower<- qu[1]-1.5*(qu[3]-qu[1])
cat("upper", upper); cat("/", "lower", lower)

## upper 33.86

## / lower 4.363

# upper inner fence
min(max(mtcars$mpg), upper)

## [1] 33.86

# lower inner fence
max(min(mtcars$mpg), lower)

## [1] 10.4

畫盒型圖 2.15 對照上述的計算結果：

boxplot(mtcars$mpg, ylim=c(0,40), yaxt='n')
axis(2, at=c(1:40, by=5), labels=c(1:40,by=5))

Figure 2.15: 盒型圖二

上述的資料並沒有極端值，我們可以觀察美國各州的面積，圖 2.16 顯示有些州的面積接近60萬平方英里，有的不到10萬平方英里：

boxplot(state.area, ylab="Area of State")

Figure 2.16: 盒型圖三

然後請輸入下列指令：

y<-state.abb

identify(rep(1, length(y)), y, labels=seq_along(y)))

然後用滑鼠對著圖形點擊，如果出現已經找到最近的點，可以按旁邊的finish，然後會呈現哪幾個觀察值被標示在圖形上面。例如第二個州是阿拉斯加，面積最大。

我們可以對類別變數繪製另一個連續變數的箱型圖，例如圖2.17：

Orange$tree <-ordered(Orange$Tree, levels=c(1,2,3,4,5))
with(Orange, boxplot(circumference ~ tree))

Figure 2.17: 盒型圖四

由上圖 2.17 可以看出，第一類型的四分位距比較小，其次是第三類型。而第四類型的中位數最大。

2.3.1 課間練習

☛請試著畫ISLR套件中Wage資料的薪水盒型圖。薪水變數是該資料的列的名稱。

---

2.4 長條圖 (Histogram)

我們可以從長條圖觀察連續變數的分佈，例如偏態以及峰值等等。 因為連續變數有一定數量的值，而非像類別變數的值是離散的，所以應該用長條圖來呈現。例如圖 2.18 顯示美國50個州的傷害罪統計，有的州落在50件以內，有的州則在100到150中間：

par(mfrow=c(1,2))
hist(USArrests$Assault, col="tomato2", main="breaks_default")
hist(USArrests$Assault, breaks = 15, col="tomato3",main="breaks_15")

Figure 2.18: 長條圖一

上圖 2.18 的指令中，我們可以指定或者不指定長條的數目。不同的長條數目呈現不同的分佈。長條數目越多可能越接近資料，但是也會帶來許多雜訊。

連續變數的每一個值的相對比例可以幫助我們了解分佈型態，因此，freq=F參數強制長條圖呈現相對的比例，而非次數。圖 2.19 顯示傷害罪的機率密度。

hist(USArrests$Assault, col="tomato4", freq = F, 
     xlab="Assault", main="Assault in 50 States", breaks = 10)

Figure 2.19: 長條圖二

R會自動挑選適合該變數分佈離散程度的寬度，breaks參數越大，寬度越小，有可能出現某一間隔沒有任何觀察值之狀況。但是寬度越大，變數的分佈越粗略。

我們用模擬的資料來表現長條圖的參數break的用途。我們從常態分佈的x抽出100個，計算平均值後記錄為y，然後重複1000次，也就是有1000個x的平均值。在圖 2.20，左邊的圖最多有10個等分，右邊的圖則是可以到50個，所以右邊的圖顯示比較多的直條。

par(mfrow=c(1,2))

y <- c()
for (i in 1:1000)
{x=rnorm(1000,0,1)
 x.sample <- sample(x, 100)
  y[i]=mean(x.sample)}

hist(y, 10, probability = T)
rug(jitter(y))
hist(y, 50, probability = T)
rug(jitter(y))

Figure 2.20: 兩個長條圖並列

在長條圖上面可以加上特定分布的曲線，例如常態分佈曲線或者是一致分佈曲線如圖 2.21：

par(mfrow=c(1,1))
set.seed(02138)
y <- c()
for (i in 1:1000)
{x=rnorm(1000,0,1)
 x.sample <- sample(x, 100)
  y[i]=mean(x.sample)}

hist(y, 100, probability = T, col="gray90")

curve(dunif(x, min=min(y), max=max(y)), add=T, col="blue", lwd=2)

curve(dnorm(x, mean=mean(y), sd=sd(y)), add=T, col="red", lwd=2)

Figure 2.21: 加上機率密度曲線的長條圖

從上面的圖可以看出模擬的資料近似常態分佈，而常態分佈曲線比一致分佈曲線更接近資料分佈。

我們再用ggplot2::diamonds中的price為例，用三種方法加上分布的曲線。由於這筆資料的單位比較大，所以先除以1000。第一種是直接加上機率密度曲線（粉紅）。第二種是加上右偏的卡方分佈曲線（藍）。第三種則是加上指數分佈曲線（黑）。指數分佈的函數為：

\[f(x|\lambda)=\lambda e^{-\lambda x}\quad x\ge 0\]

圖2.22顯示價格的機率密度：

Y<-ggplot2::diamonds$price/1000

hist(Y, freq =  F , breaks=25)
lines(density(Y), lty=2, 
      lwd=2, col="#e122aa")


curve(dchisq(x, df=4), add=T, 
      lwd=2, col="darkblue")

curve(dexp(x, rate=.25), lwd=2, add=T)

Figure 2.22: 鑽石價格機率密度

有關各種分布的語法可參考Statistics Globe這個網站。在後續的課程中我們也會介紹更多的分佈。

---

2.5 折線圖 (Line Chart)

之前在介紹散佈圖時，曾經設定參數type='n'，使觀察值不顯示出來，而除了type='n'，還有b, c, h, o, p等5個選項。而且，可以用單一變數繪圖表示分佈型態。

例如圖 2.23 顯示風速的趨勢：

with(airquality, plot(Wind, type='b'))

Figure 2.23: 空氣品質折線圖一

或者是圖 2.24 設定不同的顏色以及符號，顯示臭氧層的變化情形：

with(airquality, plot(Ozone, type='o', pch=16, 
      cex=1.2, lty=2, lwd=2, col='red'))

Figure 2.24: 空氣品質折線圖二

2.5.1 時間序列資料

資料有時間序列的型態，適用以折線圖顯示，例如LakeHuron這筆資料是時間序列資料，裡面的參數設定為type='o'。圖2.25顯示湖水的折線圖：

plot(LakeHuron, type = "o", pch=16, cex=1.2, lty=2)  ## Index plot

Figure 2.25: Lake Huron折線圖

那麼我們要如何建立時間序列資料？ts()函式可以轉換一個資料框為時間序列資料，要注意設定起始的時間點，可以以月份來切分，如果頻率為1，就是以1年為單位，如果頻率為2，就是半年，頻率為4，就是3個月：

CSV<-here::here("data","Tondutrend.csv")
trend<-read.csv(CSV, header=T, sep=",")
tonduts<-ts(trend, start=1992.6, frequency=2)
tonduts

## Time Series:
## Start = 1992.6 
## End = 2019.6 
## Frequency = 2 
##      統一 維持現狀 獨立  拒答
## 1993  2.4     66.1  7.7 23.80
## 1993 18.6     42.2  4.0 35.20
## 1994 20.4     51.4 11.4 16.70
## 1994   NA       NA   NA    NA
## 1995   NA       NA   NA    NA
## 1995   NA       NA   NA    NA
## 1996 19.1     46.1 11.1 23.80
## 1996 22.3     41.1 12.2 24.40
## 1997 21.0     48.5 14.1 16.40
## 1997 27.8     47.1 16.1  9.00
## 1998 22.8     50.8 17.9  9.00
## 1998 20.5     49.9 17.4 12.20
## 1999 20.3     50.5 20.0  9.20
## 1999 18.8     53.8 18.9  8.50
## 2000 13.3     58.8 17.8 10.10
## 2000 20.8     52.5 16.2 10.60
## 2001 26.8     51.1  9.0 13.00
## 2001 23.8     55.5 13.3  7.40
## 2002 26.2     52.5 15.4  6.00
## 2002 21.5     53.5 17.9  7.10
## 2003 21.6     52.9 18.8  6.80
## 2003 17.6     56.6 19.2  6.60
## 2004 14.9     55.9 21.6  7.70
## 2004 12.8     60.2 19.6  7.40
## 2005 14.1     57.7 20.7  7.50
## 2005 15.7     59.3 19.8  5.20
## 2006 16.2     58.4 20.1  5.20
## 2006 14.9     60.9 18.5  5.70
## 2007 15.9     60.7 18.6  4.80
## 2007 18.2     53.3 23.1  5.40
## 2008 12.9     58.2 20.8  8.00
## 2008 11.5     59.6 20.4  8.50
## 2009 10.5     59.6 25.7  4.10
## 2009 10.2     63.1 21.3  5.50
## 2010 10.6     63.9 19.6  5.90
## 2010 12.3     61.0 22.8  3.80
## 2011 10.0     64.6 21.0  4.40
## 2011 11.1     61.3 23.0  4.60
## 2012 10.8     63.7 19.1  6.40
## 2012 10.7     64.8 19.1  5.40
## 2013 11.1     62.8 19.6  6.50
## 2013 12.3     59.2 23.3  5.20
## 2014 11.6     61.3 21.7  5.40
## 2014 10.9     60.8 22.9  5.40
## 2015  8.8     61.6 23.6  5.90
## 2015  9.6     61.6 20.5  8.40
## 2016 10.8     60.6 22.2  6.40
## 2016 10.0     61.3 22.8  5.90
## 2017 11.4     60.1 22.5  5.90
## 2017 13.2     58.7 23.1  5.00
## 2018 14.4     59.8 19.9  5.90
## 2018 16.5     58.1 20.0  5.50
## 2019 17.5     57.9 19.5  5.20
## 2019 10.6     57.6 25.8  6.00
## 2020  7.5     57.6 28.0  6.98

相對於其他軟體，R的優點在於允許時間點有缺漏，畫出來的圖形就有缺少某個或是某些時間點，不會強制連接前後兩個時間點。接下來是畫折線圖。圖 2.26 顯示2004年以來統獨趨勢一直相當穩定，但是在2016年之後偏統的支持度似乎上升，到了2019年又下降，而偏獨立的民意上升：

# Margin
par(xpd=NA, mar=par()$mar+c(2.5, 0, 0, 0),family='HanWangKaiMediumChuIn')
# Plot
plot(tonduts, plot.type=c("single"), lty=c(1,2,3,2),ylab="%",xlab=NULL,pch='1', lwd=3,frame.plot=F,col=c("gray20","gray60", "black", "gray80"),xaxt="n", main="台灣民眾的統獨立場, 1992.6-2020.12")
axis(1, at=seq(1992,2020,by=2))
axis(2, at=seq(10,70,by=10))
legend("bottomright", c("統一","維持現狀"),inset=c(0.35, -0.6), col=c("gray20","gray60"),lty=c(1,2), bty='n', lwd=3)
legend("bottomright", c("獨立","無反應"), inset=c(0, -0.6), col=c("black", "gray80"),lty=c(3,2), bty='n', lwd=3)
text(2000, 30, paste("第一次政黨輪替"))

Figure 2.26: 統獨趨勢圖

在圖形 2.26中，我們設定xpd=NA，允許圖形超過界線，而且設定圖形的邊線可以到原來區域的底部2.5個文字。這是因為折線有四條，需要一定的空間容納圖例，而為了避免與圖形重疊起見，我們把圖例放在X軸右下方，所以除了要設定參數"bottomright"，還要設定inset(0.35, -0.6)以及inset(0, -0.6)。 除了設定圖例出現的位置，還可以設定bty="n"，強制圖例沒有外框。

有時候R會出現figure margin too large的錯誤訊息，這時候請調整一下par(mar=c())，例如par(mar=c(-0.2,-0.2,0,0))

最後，我們加上一點註解文字到圖形上面，當然，我們可以加上線條或是箭頭，請大家研究一下segment這個指令。相信這個折線圖會讓你的報告比其他軟體畫的圖看起來更專業。

這張圖裡面有中文，所以練習時，我們需要設定中文字型，不然會顯示亂碼。如果你打開字體簿，可以找PostScript名稱，然後在par()裡面設定參數family=""，就像底下的語法。

library(here);library(foriegn)
CSV<-here("data","Tondutrend.csv")
trend<-read.csv(CSV, header=T, sep=",")
tonduts<-ts(trend, start=1992.6, frequency=2)
par(xpd=NA, mar=par()$mar+c(2.5, 0, 0, 0), 
    family='HanWangKaiMediumChuIn')

另一個做法是找字體的檔案，設定字體，就像底下的語法。

library(here);library(foriegn);library(showtext)
showtext.auto()
font.add("細明體", "MingLiu.ttf")
par(xpd=NA, 
    mar=par()$mar+c(0, 0, 0, 0), 
    family='細明體')

2.5.2 課間練習

☛請畫圖表示MASS::accdeaths的趨勢。

---

2.6 特殊點狀圖

symbols()可以產生指定形狀、大小的散佈圖，而且可以進一步根據另一個變數的觀察值，調整散佈點的大小，例如我們想顯示兩個變數之間關係，加上另一個連續變數的大小程度，見散佈圖2.27：

with(anscombe, symbols(x2, y1, circles=anscombe$x1,   inches=0.2, fg='blue'))

Figure 2.27: 特殊點狀圖一

可以看出有些圈圈比較大，有些比較小。再用美國的犯罪統計舉例，見圖2.28：

with(USArrests, symbols(Murder, Assault, circles=UrbanPop, inches=0.12, bg="red"))

Figure 2.28: 特殊點狀圖二

上圖顯示，雖然搶案發生率與暴力犯罪率成正比，都市化人口越多的地方，也可能有更多的類似案件。我們可以進一步用迴歸統計確認。

m.arrest<-with(USArrests, lm(Assault ~ Murder+UrbanPop))
summary(m.arrest)

## 
## Call:
## lm(formula = Assault ~ Murder + UrbanPop)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -107.78  -25.99   -1.97   22.49  111.82 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -23.620     33.217   -0.71    0.481    
## Murder        15.071      1.572    9.59  1.2e-12 ***
## UrbanPop       1.175      0.473    2.48    0.017 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 47.8 on 47 degrees of freedom
## Multiple R-squared:  0.684,  Adjusted R-squared:  0.671 
## F-statistic:   51 on 2 and 47 DF,  p-value: 1.69e-12

2.6.1 課間練習

☛請畫圖表示MASS::Boston這筆資料中，與生師比(ptratio)與房價中位數(medv)的關係，並且考慮低社會地位人口比例(lstat)的作用。

2.7 輸出圖形

• 我們可以輸出圖形到指定的路徑，方便放在報告當中。例如我們產生一個圓餅圖，然後輸出為png檔，在執行繪圖指令前，要先執行輸出圖形的指令，例如：

#設定路徑
setwd(here::here())
setwd(here::here('Fig'))
#設定圖形輸出的檔案名稱與大小
png("pie.png", width=4, height = 6, units = 'in', res=300)
#繪圖
pie(table(ISLR::Auto$origin))
#結束
dev.off()

quartz_off_screen 2

---

3 圖形進階(ggplot2)

接下來我們介紹用ggplot2繪圖。基本的指令為：

ggplot(data, aes(x, y, group,...))+ geom_object()+theme()

• geom_object可以是geom_bar（直方圖）、geom_point（散佈圖）、geom_line（折線圖）等等。
  
• 繪圖的對象可以是資料框或是已經轉為表格的現成資料。
  
• 繪圖之前應該特別注意變數的型態。
  
• 這個網頁有許多ggplot2的例子，有興趣的同學可參考。

---

3.1 直方圖 (Bar plot)

變數屬於數字、類別或是字串都可以以直方圖表示。例如圖 3.1顯示單一變數的直方圖：

ggplot(ChickWeight, aes(x=Diet), stat='bin') +
  geom_bar(fill='red', col='blue') +
  ggtitle('Diet')

Figure 3.1: ggplot2直方圖一

或者是先轉換成表格，但是要設stat="identity"，畫成3.2：

stu<-here::here('data','studentsfull.txt')
students<-read.table(stu, header=T, sep='')
stu.t <- as.data.frame(table(students$Department, dnn=c("Dep")))
stu.t

##          Dep Freq
## 1  Aerospace    4
## 2  Chemistry    3
## 3  Economics    4
## 4    English    3
## 5 Journalism    4
## 6  Mechanics    5
## 7    Physics    3

g1 <- ggplot(stu.t, aes(x=Dep, y=Freq)) + 
geom_bar(stat = 'identity', fill="#4411FF33") + ylab('Count')
g1

Figure 3.2: ggplot2直方圖二

我們也可以直接用現成的次數分配畫圖。先輸入四家公司的市值：

dt <- scan(what=list(company="character", 
                     marketvalue="numeric"))
#Apple 851
#Microsoft 703
#Amazon 701
#Facebook 464
#DT<-data.table::setDT(dt)

然後我們轉成資料框，再畫圖 3.3。我們可以用aes給每一個直方一個顏色，fill=company。aes的功能類似相對應的函數，例如company是一個變數，而fill=company就是根據這個變數選擇顏色。

dt <- data.frame(company=c("Apple", "Microsoft",
                           "Amazon", "Facebook"),
                 marketvalue=c(851,703,701,464))


DT <- dt %>% mutate(company=as.factor(company),
                    marketvalue=as.numeric(marketvalue))
ggplot(DT, aes(x=reorder(company, -marketvalue), 
         y=marketvalue)) +
  geom_bar(aes(fill=company),stat="identity") +
       theme_classic()

Figure 3.3: 次數分配轉成直方圖

在畫圖 3.3 時要注意，X軸需要「類別」，而Y則是數值。 我們在ggplot2裡面用reorder()函數，改變X軸的類別順序為由大排到小。

3.1.1 課間練習

☛請由小而大，呈現reshape2的tips資料裡面，週四到週日的每天(day)收到的小費(tip)的平均數。請問哪一天收到小費的平均數最大？

3.1.2 加上百分比到Y軸

我們可以用aes(y=(..count..)/sum(..count..))轉化類別為百分比，如圖 3.4。

ggplot(data=tips, aes(x=time)) +
     geom_bar(aes(y=(..count..)/sum(..count..)),
    stat="count", fill="#ee0011") +
      scale_y_continuous(label = scales::percent) 

Figure 3.4: 加上百分比的直方圖

3.1.3 加上數據到直方圖

用已經統計好的數字當作標記(label)，而該數字加上或是減去某個單位當作y，用 geom_text() 標上數據，例如圖 3.5 顯示各系的人數：

g1 + geom_text(data=stu.t, aes(x=Dep, label=Freq, y=Freq+0.2), size=8)

Figure 3.5: 加上資料標籤的直方圖

3.1.4 應用stat_count()

stat_count()可以跟geom_bar()互換，一樣產生直方圖，例如圖 3.6 為加上類別的個案數的直方圖：

salary<-alr4::salary
ggplot(salary, aes(x=rank)) + 
  stat_count(geom='bar', fill='#0aaee1') +
  geom_text(aes(label=..count..), stat='count', colour='white', vjust=1.6)

Figure 3.6: ggplot2直方圖三

3.1.5 兩個類別變數的直方圖

我們想用直方圖表現兩個變數之間的關係，例如在學生成績這個資料中，我想想知道每一個系考試超過70分的比例，也就是系與學生是否通過兩個變數之間的關係。

首先，我們用 dplyr::mutate 產生一個新的變數 Pass：

full<-here::here('data','studentsfull.txt')
students<-read.table(full, header=T, sep='')
stu <- dplyr::mutate(students, Pass=(Score>70))
kableExtra::kable_styling(knitr::kable(stu))

ID	Name	Department	Score	Gender	Pass
10322011	Ariel	Aerospace	78	F	TRUE
10325023	Becky	Physics	86	F	TRUE
10430101	Carl	Journalism	69	M	FALSE
10401032	Dimitri	English	83	M	TRUE
10307120	Enrique	Chemistry	80	M	TRUE
10207005	Fernando	Chemistry	66	M	FALSE
10305019	George	Mechanics	75	F	TRUE
10305022	Howell	Mechanics	81	M	TRUE
10305029	Ian	Mechanics	60	M	FALSE
10305031	Julio	Mechanics	89	M	TRUE
10322014	Kaori	Aerospace	82	F	TRUE
10425026	Luke	Physics	88	M	TRUE
10401022	Miguel	English	92	M	TRUE
10501006	Neo	English	77	M	TRUE
10321010	Olivia	Economics	85	F	TRUE
10321011	Peter	Economics	88	M	TRUE
10405017	Qing	Mechanics	88	F	TRUE
10422007	Ricky	Aerospace	91	M	TRUE
10422008	Seiko	Aerospace	80	F	TRUE
10430005	Terresa	Journalism	62	F	FALSE
10530009	Usla	Journalism	87	F	TRUE
10421001	Vivian	Economics	70	F	FALSE
10307018	Wendy	Chemistry	85	F	TRUE
10425003	Xing	Physics	93	M	TRUE
10221030	Yoko	Economics	66	F	FALSE
10430015	Zoe	Journalism	92	F	TRUE

然後就可以畫成並列的直方圖區分通過與不通過，如圖 3.7：

ggplot(stu, aes(x=Department, fill=Pass)) + 
          geom_bar(position='dodge')

Figure 3.7: 並列直方圖

如果要讓色塊表示條件機率，我們可以用 summarize 產生 Count這個變數，代表以系為單位，每個系通過以及沒通過的個案數。然後用 mutate 產生以系所為邊際機率的通過與不通過的條件機率。

library(dplyr)
stu <- students %>%  mutate(Pass=(Score>70))
stu.ag <- summarize(group_by(stu, Department, Pass), Count=n())
stu.ag <- stu.ag %>% mutate(Pct=Count/sum(Count))
kableExtra::kable_styling(knitr::kable(stu.ag))

Department	Pass	Count	Pct
Aerospace	TRUE	4	1.0000
Chemistry	FALSE	1	0.3333
Chemistry	TRUE	2	0.6667
Economics	FALSE	2	0.5000
Economics	TRUE	2	0.5000
English	TRUE	3	1.0000
Journalism	FALSE	2	0.5000
Journalism	TRUE	2	0.5000
Mechanics	FALSE	1	0.2000
Mechanics	TRUE	4	0.8000
Physics	TRUE	3	1.0000

比較表@ref(tab:students.summarize) 跟表@ref(tab:students.mutate)，可以發現前者是後者的統計，也就是統計有幾個人通過，例如航太系有4人通過，化學系有2人通過，1人沒通過等等。

畫成圖表示上列數據如圖 3.8：

ggplot(stu.ag, aes(x=Department, y=Pct, fill=Pass)) +
   geom_bar(stat='identity') +
  scale_y_continuous(label = scales::percent)

Figure 3.8: ggbarplot2堆疊直方圖

• 要注意Y軸是Pass的各類別所佔的比例，因此fill=Pass

如果只想顯示通過的比率，可篩選掉未通過的資料，然後畫成圖 3.9：

newstu<-stu.ag[stu.ag$Pass==TRUE,]
ggplot(newstu, aes(x=Department, y=Pct, fill=Pass)) +
   geom_bar(stat='identity') +
  scale_y_continuous(label = scales::percent)

Figure 3.9: ggplot2堆疊直方圖

3.1.6 課間練習

☛請用ggplot2::diamonds這筆資料，畫出切割(cut)以及成色(color)兩個變數的堆疊直方圖，呈現不同的切割品質的鑽石中，成色所佔的比例。

---

如果我們只想知道某一個類別變數在另一個類別變數的各類別中的細分類，qplot可以很快地產生圖形。以鑽石的切割以及成色為例，圖 3.10 顯示每一個切割等級的個數以及其中各種成色的數目：

qplot(cut, data = diamonds, 
      geom = "bar",  fill = color, stat="count")

Figure 3.10: qplot直方圖

循著之前的方法，以N=n()計算每個成色在每個類別中的原始數目，然後把N放在Y軸。請見圖 3.11。

df<-ggplot2::diamonds
df <- df %>% select(cut, color)
newdf <-reshape2::melt(df, id.vars=c("cut"))
colnames(newdf)<-c("cut", "color.g", "color")

qdf <- summarize(group_by(newdf, cut, color), N=n())  
g1 <- ggplot(qdf, aes(x=cut, y=N, fill=color)) +
           geom_bar(stat="identity")
g1 + scale_fill_grey(start = 0.2, end = 0.9) +
     theme_bw()

Figure 3.11: 灰階原始數目直方圖

---

3.2 長條圖 (Histogram)

ggplot2的長條圖使用的函數為geom_histogram()跟直方圖類似，但是適用在連續變數。參數stat='bin'指定每一長條的寬度，寬度越寬，越容易黏合在一起，但是也可能越不容易分辨分佈的型態，例如圖 3.12 顯示mtcars資料中，mpg變數的長條圖：

ggplot(mtcars, aes(x=mpg)) + 
 geom_histogram(stat='bin', binwidth=1, fill="#a310ec")

Figure 3.12: ggplot2長條圖一

• binwidth 越接近1，長條的寬度越大，長條的數目也越少，也就越集中。所以要挑選適當的寬度，才能正確地表現出資料的分佈。

長條圖可以轉換成密度圖，而且可以按照對應的類別變數並列在同一張圖，以圖 3.13為例，geom_density()可以顯示機率密度，並且同時顯示三個類別：

ggplot(mtcars, aes(x=wt, fill=as.factor(am))) +
  geom_density(position="identity", alpha=.4)

Figure 3.13: ggplot2長條圖二

3.2.1 課間練習

☛請畫出UsingR::batting這筆資料中的「長打率」分佈圖。長打率指的是：

\[\frac{H+2\times Double+3\times Triple+4\times HR}{AB}\]

---

3.3 散佈圖 (scatter plot)

如果需要表現兩個變數之間的關係，散佈圖可以顯示兩個變數對應的位置。散佈圖用的函數是geom_point()。例如圖3.14顯示車的重量與mpg的關係：

ggplot(mtcars, aes(x=wt, y=mpg)) +
    geom_point()

Figure 3.14: ggplot2散佈圖一

我們可以控制散佈圖的符號的形狀、大小、顏色等等。例如圖 3.15：

set.seed(02138)
ggplot(data=data.frame(x=rnorm(100, 1, 1.5),y=rnorm(100,0,1)),
       aes(x=x,y=y))+
  geom_point(col="#b10a2c", shape=5, size=2)

Figure 3.15: ggplot2散佈圖二

3.3.1 加上迴歸線

散佈圖可以加上迴歸線表現兩個變數的相關，如圖3.16：

ggplot(mtcars, aes(x=wt, y=mpg)) +
    geom_point(col='red') +
    geom_smooth(method="lm", se=F, col='blue')

Figure 3.16: ggplot2散佈圖加上迴歸線

也可以用非線性迴歸線表示兩者的相關，如圖3.17：

ggplot(Orange, aes(x=age, y=circumference)) +
   geom_point(aes(col=Tree), size=3) +
   geom_smooth(method="loess", se=F) +
   theme_bw()

Figure 3.17: 散佈圖加上無母數迴歸線

3.3.2 課間練習

☛請用散佈圖加迴歸線表示長打率跟三振率(SO/AB)的關係：

---

3.4 三個變數之散佈圖

兩個變數的散佈圖可表示兩個變數之間的相關，例如圖 3.18：

ggplot(anscombe, aes(x=x1, y=y1)) +
  geom_point(size=3) 

Figure 3.18: 三個變數散佈圖一

如果加上另一個變數，可以表示當另一個變數等於某一個類別時，兩個變數是否仍然維持一樣的關係？我們用shape設定形狀隨著類別變數而改變，然後隨著類別變數改變顏色如圖 3.19：

n <- nrow(anscombe)
anscombe$r <- rep(0, n)
anscombe$r[anscombe$x1>8]<-1
sc1<-ggplot(anscombe, aes(x=x1, y=y1, shape=factor(r))) 
sc1 + geom_point(aes(color=factor(r)) , size=3) 

Figure 3.19: 三個變數散佈圖二

也可以隨著連續變數而改變顏色深淺如圖 3.20：

sc2<-ggplot(anscombe, aes(x=x1, y=y1)) 
sc2 + geom_point(aes(color=x2) , size=3) 

Figure 3.20: 三個變數散佈圖三

3.4.1 課間練習

☛reshape2裡面的 french_fries 針對油炸與薯條的味道進行實驗。在這筆資料裡有幾個測量薯條口味的變數，請找兩個變數，畫出散佈圖，但是同時要顯示三個實驗組別的分佈。

---

3.5 盒型圖 (Box plot)

盒型圖可以表現一個變數或者是一個類別變數對應另一個連續變數的分佈，不過ggplot2無法直接畫出一個連續變數的盒型圖，所以我們先創造一個只有一個值的變數，然後對應我們要顯示的連續變數，例如圖 3.21：

mtcars <- mutate(mtcars, X=1)
ggplot(mtcars, aes(x=X, y=mpg)) +
     geom_boxplot() +
  labs(x="",y='mpg') +
  stat_summary(fun.y=median, geom="point", shape=16, size=2) +
  theme_bw()

Figure 3.21: 盒型圖一

盒型圖可以用來比較一個類別變數對應的連續變數的分佈，例如我們想要顯示不同的氣缸數（cyl）的馬力（hp），可以把y設定為（hp），然後把x, group設定為（cyl），如圖 3.22：

g1 <-ggplot(mtcars, aes(x=cyl, y=hp, group=cyl)) +
  geom_boxplot()  +
  labs(x='Cylinder', y="Horse Power") +
  stat_summary(fun.y=median, geom="point", shape=16, size=2) +
  theme_bw()
g1

Figure 3.22: 盒型圖二

由於 cyl 是數字變數，我們發現有一些我們不需要的類別，沒有資料與其對應，因此我們用 scale_x_discrete加以調整，也就是只留下需要的類別數字，然後加上所需要的標籤如圖 3.23：

g1 + scale_x_discrete(limit = c(4, 6, 8),
       labels = c("4","6","8")) + theme_bw()

Figure 3.23: 盒型圖三

類似的技巧可以適用在Y軸，例如我們只想顯示50到200的馬力，超過200以上的資料就會被去掉如圖 3.24：

g1 + scale_y_continuous(limits = c(50,200)) +
       scale_x_discrete(limit = c(4, 6, 8),
              labels = c("4","6","8")) +
              theme_bw()

Figure 3.24: 盒型圖四

3.5.1 課間練習

☛請分別畫出 UsingR::movies} 這筆資料中，current、previous、gross這三個變數的盒型圖以及只有current、previous這兩個變數的盒型圖。

3.6 甜甜圈圖

過去有所有的圓餅圖，甜甜圈圖可以顯示資料的次數分佈，例如我們想知道ggplot2裡面的的diamonds資料中的cut的次數分配，先針對這個變數計算個數：

dt <-group_by(diamonds, cut)
dat = summarize(dt, count=n())
kableExtra::kable_styling(knitr::kable(dat))

cut	count
Fair	1610
Good	4906
Very Good	12082
Premium	13791
Ideal	21551

再計算比例、排序以及計算累積百分比：

# Add addition columns, needed for drawing with geom_rect.
dat$fraction = dat$count / sum(dat$count)
dat = dat[order(dat$fraction), ]
dat$ymax = cumsum(dat$fraction)
dat$ymin = c(0, head(dat$ymax, n=-1))
kableExtra::kable_styling(knitr::kable(dat))

cut	count	fraction	ymax	ymin
Fair	1610	0.0298	0.0298	0.0000
Good	4906	0.0910	0.1208	0.0298
Very Good	12082	0.2240	0.3448	0.1208
Premium	13791	0.2557	0.6005	0.3448
Ideal	21551	0.3995	1.0000	0.6005

然後畫圖 3.25：

p1 = ggplot(dat, aes(fill=cut, ymax=ymax, ymin=ymin, xmax=4, xmin=3)) +
     geom_rect() +
     coord_polar(theta="y") +
     xlim(c(0, 4)) +
     theme(panel.grid=element_blank()) +
     theme(axis.text=element_blank()) +
     theme(axis.ticks=element_blank()) +
     annotate("text", x = 0, y = 0, label = "Diamond Cut") +
     labs(title="") +
     xlab("") +
     ylab("")
p1

Figure 3.25: 甜甜圈圖一

3.6.1 課間練習

☛請嘗試把UsingR::firstchi這筆資料分成四組，然後畫出甜甜圈圖表示四組的相對大小。

---

3.7 折線圖

之前我們用R的基本指令對airquality這筆時間序列資料畫折線圖，例如圖 3.26：

ggplot(airquality, aes(x=Day, y=Wind)) +
  geom_line()

Figure 3.26: ggplot折線圖一

這個圖看起有很多重疊的地方，原因是日期是從1日到30日或是31日循環出現，所以我們可以創造一個新的變數，可以顯示時間序列，例如圖 3.27：

N <- summarize(airquality, n())
airquality <- mutate(airquality, index=seq_along(1:N[,1]))
ggplot(airquality, aes(x=index, y=Wind)) +
  geom_line() 

Figure 3.27: ggplot折線圖二

我們示範如何把原來的月份與日變數轉換成日期：

airquality$Date <- paste0(airquality$Month, "/", airquality$Day) 
airquality$Date  <- as.Date(airquality$Date, format="%m/%d")

然後用日期畫折線圖 3.28：

ggplot(airquality, aes(x=Date, y=Wind)) +
  geom_line() 

Figure 3.28: ggplot2折線圖三

進一步可以顯示不只一項的折線圖，但是我們需要轉置資料為長表，在轉置之前先選擇兩個變數加上日期：

airquality.n <- select(airquality, Date, Wind, Temp)
library(reshape2)
airquality.l<-melt(airquality.n, id.vars=c('Date'))
head(airquality.l)

##         Date variable value
## 1 2021-05-01     Wind   7.4
## 2 2021-05-02     Wind   8.0
## 3 2021-05-03     Wind  12.6
## 4 2021-05-04     Wind  11.5
## 5 2021-05-05     Wind  14.3
## 6 2021-05-06     Wind  14.9

在轉換之後，出現一個新變數value，這個新變數包含了Wind, temp兩個變數的值，這樣我們才能一次畫兩個變數的值。

由於這兩個變數的大小範圍接近，可以畫在一張圖上面。例如圖 3.29：

ggplot(airquality.l, aes(x=Date)) +
  geom_line(aes(y=value, col=variable)) 

Figure 3.29: ggplot2折線圖四

3.7.1 課間練習

☛請畫出alr4::ftcollinssnow這筆資料中1到6月下雪(Late)量的折線圖。

3.8 特殊圓點圖

ggplot2可以應用在比較迴歸係數的大小以及標準誤的範圍，方便讀者瞭解模型的估計結果。例如我們估計有三個自變數的迴歸模型：

m1 <- with(mtcars, lm(mpg ~ wt + hp + am))
summary(m1)

## 
## Call:
## lm(formula = mpg ~ wt + hp + am)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.422 -1.792 -0.379  1.225  5.532 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 34.00288    2.64266   12.87  2.8e-13 ***
## wt          -2.87858    0.90497   -3.18  0.00357 ** 
## hp          -0.03748    0.00961   -3.90  0.00055 ***
## am           2.08371    1.37642    1.51  0.14127    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.54 on 28 degrees of freedom
## Multiple R-squared:  0.84,   Adjusted R-squared:  0.823 
## F-statistic:   49 on 3 and 28 DF,  p-value: 2.91e-11

然後我們建立一個資料框，裡面包含係數以及對應的標準誤，但是為了分析方便，我們去掉截距。

同時，我們產生兩個變數，第一個是標準誤的下限，第二是標準誤的上限：

df <- data.frame(coef=coef(m1)[-1], se=coef(summary(m1))[-1, 
                                    "Std. Error"] )

dt <- mutate(df, lower=coef-qt(0.975, 
            nrow(mtcars))*se, upper=coef+qt(0.975, 
                        nrow(mtcars))*se)

一切就緒後，我們利用geom_segment這個參數，產生線段，原理是x到xend，以及y到yend分別產生線段如圖 3.30：

s1 <- ggplot() + 
  geom_segment(data=dt, 
  mapping=aes(x=row.names(dt),  y=lower, 
              xend=row.names(dt), yend=upper), 
                size=2, color="blue") 
s1

Figure 3.30: ggplot圓點圖一

然後加上代表係數的圓點，並加上係數的名稱，以及把X軸與Y軸對調，就大功告成。圖 3.31 顯示是否自動排檔對於依變數的影響相對來得大，但是有可能等於0，因為標準誤的下限低於0。

s1 + geom_point(data=dt, 
                aes(x=row.names(dt), y=coef), size=4, shape=21, fill="white") +
    scale_x_discrete(breaks=c('wt','hp','am'),
        labels=c("Weight", "Horse Power", "Auto")) +
    xlab("") +
    ylab("Estimates")  +
    coord_flip()  +  
 theme(axis.text.y = element_text(face="bold", color="#993333",  size=10))

Figure 3.31: ggplot圓點圖二

這個係數圖雖然已經相當美觀，但是仔細觀察，會發現係數的順序反過來了。要改變這個情形，可以用rev()這個函數改變向量的元素順序，例如：

V<-c(50, 100, 30, 200)
rev(V)

## [1] 200  30 100  50

3.8.1 課間練習

☛請把係數名稱跟係數順序根據 3.30 得到的順序排列一遍：

tmp <- dt[order(rev(row.names(dt))),]
s2=ggplot() + 
    geom_point(data=tmp, aes(x=row.names(tmp), 
                y=rev(coef)), size=4, shape=21, fill="#e21b33") +
    scale_x_discrete(breaks=c("wt","hp","am"),
        labels=c( "Auto", "Horse Power","Weight")) +
    ylab("") +
    xlab("")  +
    coord_flip()  +  
    theme(axis.text.y = element_text(face="bold", color="#00ABFE", size=14) )
s2+geom_segment(data=dt, mapping=aes(x=row.names(dt),  y=rev(lower), xend=row.names(dt), yend=rev(upper)), 
                size=2, color="#FABEC1") 

Figure 3.32: ggplot2係數圖

3.9 特殊風格

除了上述的圖形，ggplot2還有更多有趣的圖形，也允許使用者設定圖形風格、資料點顏色、副標題、圖例的標題等等，請同學上網多多參考他人的例子。例如圖3.33顯示經濟學人雜誌風格的直方圖：

library(ggthemes)
# Number of Cases
dt <- data.frame(Year=c(2013, 2014, 2015, 2016, 2017),
                 Amount=c(65, 56, 78, 85, 79))
# Economist theme
p1<-ggplot() + theme_economist() + scale_fill_economist() +
geom_bar(aes(x=Year, y=Amount),
         data=dt, stat='identity')
#data label
p3 <- p1+ theme(axis.title= element_text(color="blue", size=14, face="bold"),
        axis.text = element_text(size=14),
        legend.text=element_text(size=18))+
   ggtitle("Number of Cases") +
      labs(y="N", subtitle="Children, Spouses, Financial, etc.")
p3

Figure 3.33: 主題風繪圖

3.10 結合統計與圖形

• 之前我們的圖形顯示資料的原始測量，但是我們可以直接呈現統計值，例如用直方圖 3.34 比較不同切割等級的鑽石的平均價格：

ggplot(diamonds, aes(cut, price)) +
       geom_bar(stat = "summary_bin", fun.y = "mean", aes(fill=cut))

Figure 3.34: 平均值直方圖

4 lattice 繪圖

lattice套件用格子(trellis)圖形來產生高品質的資料視覺化，尤其是用在多變量的資料。

Deepayan Sakar在2008年出版「Lattice: Multivariate Data Visualisation in R」，詳細地介紹lattice套件的功能。有興趣的同學可以閱讀。

我們很快地介紹lattice套件的直方圖、長條圖、雙變數直方圖、散佈圖、折線圖、雙變數長條圖、箱型圖、點狀圖功能等等。

4.1 直方圖 (Bar chart)

我們用barchart()畫直方圖顯示星際大戰中的演員膚色如圖 4.1：

library(lattice)
barchart(starwars$skin_color,col="gray60",
       scale=list(cex=0.5, relation="free"),
         aspect=.8)

Figure 4.1: lattice直方圖一

lattice的直方圖是橫的，適合類別多的變數，例如我國面對最重要的問題，或者最常看的網站，如果用橫的直方圖會比較容易呈現。

4.1.1 課間練習

☛請嘗試對reshape2::mpg這筆資料的class做一個橫的直方圖，並且指出哪一種類型的車最多。

4.2 長條圖 (histogram)

長條圖可用來描述中央趨勢，例如圖 4.2顯示噴發時間的分佈：

histogram(faithful$eruptions, col='#a01eb2')

Figure 4.2: lattice長條圖一

histogram跟R基礎指令hist函數得到的圖形相比，多了一個外框。

4.3 雙變數直方圖 (barchart)

lattice並沒有單變數的直方圖，以長條圖代替直方圖。但是比較類別變數與連續變數之間的關係時，可以用直方圖，並且可以控制另一個類別變數。例如圖 4.3顯示品種跟顏色的關係：

library(MASS)
data(crabs)
crabs$sp<-factor(crabs$sp, labels=c("Blue",
                                     "Orange"))
crabs$sex<-factor(crabs$sex, labels=c("Female",
                                     "Male"))
barchart(sp ~ BD|sex, data=crabs, 
         col=RColorBrewer::brewer.pal(n=4, 'RdGy'), 
            main=paste0("Crab's Body Depth(mm)"))

Figure 4.3: lattice直方圖二

4.4 表格資料點狀圖 (Dotplot)

我們可以直接對一個表格資料畫圖，表示兩個變數之間的關係。例如使用VADeaths這筆表格資料，畫圖如圖 4.4：

lattice::dotplot(VADeaths, group=F, type="o")

Figure 4.4: lattice點狀圖一

以上的圖 4.4顯示，四群觀察值的年齡對應的死亡率不太一樣。

如果要畫在同一個座標軸上，需要轉換資料為資料框，如圖 4.5：

vad <- as.data.frame.table(VADeaths)
names(vad) <- c("age", "demographic", "deaths")
head(vad, n=3)

##     age demographic deaths
## 1 50-54  Rural Male   11.7
## 2 55-59  Rural Male   18.1
## 3 60-64  Rural Male   26.9

dpt<-dotplot(age ~ deaths, vad, group=demographic,
           type = "o")
update(dpt, auto.key=list(points=T, lines=T))

Figure 4.5: lattice點狀圖二

在上面的例子中，

• 先用as.data.frame.table()把表格轉成資料框
• 再用dotplot()畫出折線圖
• 最後加上圖例

4.5 統計資料點狀圖

我們想知道每一類型的樹的平均樹圍，並且畫圖。首先用tapply計算平均值。然後用dotplot畫圓點表示平均樹圍。先用tapply()函數產生每一群樹木的平均值，然後用這個平均值的表格的名稱當作Y軸。請見圖 4.6。有關apply, tapply, lapply, sapply等一連串的功能，請翻閱個人的著作第5章。

orange.mean <- tapply(Orange$circumference,
                    Orange$Tree,  mean)
dotplot(names(orange.mean) ~ orange.mean,
        aspect = .5, 
ylab = "Group",
xlab = "Mean Circumference")

Figure 4.6: lattice點狀圖一

運用dplyr的summarise指令，也可以計算平均值，記得要轉換數字變數的Tree為字串或者類別變數。請見圖 4.7。

o.mean<-Orange %>% group_by(Tree) %>% summarise(mean(circumference))
dotplot(as.factor(Tree) ~ `mean(circumference)`, aspect=0.5, data=o.mean,
ylab = "Group",
xlab = "Mean Circumference")

Figure 4.7: lattice點狀圖二

• 要注意dotplot不能用dplyr的pipeline，而是要在函數內設定data=。

我們可以用reorder把資料的順序根據平均值由小而大反過來，然後當作x，如圖 4.8：

o.mean<-Orange %>% group_by(Tree) %>% summarise(avg=mean(circumference))

dotplot(reorder(rownames(o.mean),-avg) ~ avg, aspect=0.5, data=o.mean,
ylab = "Group",
xlab = "Mean Circumference",
cex=1.5)

Figure 4.8: lattice點狀圖三

以上語法可參考這篇文章有更詳細的說明。

4.6 散佈圖 (xyplot)

兩個連續變數可以用散佈圖表示相對位置，例如圖4.9顯示E與NOx之間的關係：

with(lattice::ethanol, xyplot(NOx ~ E))

Figure 4.9: lattice散佈圖一

xyplot()用途與plot()一樣，只是內建的顏色不同。但是xyplot()可控制另一個變數，例如圖 4.10分成5個子圖：

with(lattice::ethanol, xyplot(NOx ~ E|C, col='#FF3300'))

Figure 4.10: lattice散佈圖二

4.7 類別變數與連續變數的長條圖

在lattice，類別變數對應連續變數的密度長條圖，可以一次展示在同一個圖上面。用以下指令畫圖4.11：

histogram(~temp|origin, data=nycflights13::weather, 
          type="density",   layout=c(1,3))

Figure 4.11: lattice長條圖

4.8 類別變數與連續變數的箱型圖

箱型圖也可以表示連續變數在每一個類別上面的分佈，例如圖 4.12：

lattice::bwplot(as.factor(cyl) ~ mpg, data=mtcars)

Figure 4.12: lattice箱型圖

4.9 折線圖

我們可以用xyplot畫時間趨勢的折線圖。MASS::accdeaths是時間序列的資料。圖 4.13 顯示1973年到1979年的意外死亡人數。

lattice::xyplot(MASS::accdeaths)

Figure 4.13: lattice折線圖一

圖 4.13顯示好幾年的趨勢，如果觀察MASS::accdeaths，可以發現年代在左邊、月份在上面，資料點從1973年1月排到1978年12月。也就是說，我們可以拆開6個年度的資料，然後分別畫趨勢圖。

TSstudio這個套件中，有 ts_reshape 這個函數，我們用這個函數轉換資料為長表。

tsp<-TSstudio::ts_reshape(MASS::accdeaths,
                  type = 'long')
head(tsp, n=5)

##   year month value
## 1 1973     1  9007
## 2 1973     2  8106
## 3 1973     3  8928
## 4 1973     4  9137
## 5 1973     5 10017

轉好資料後，我們用xyplot()畫圖，方法為\(y \sim x\)，因為我們有年代這個變數，當做一個中介變數(conditioning variable)，函數寫成\(y \sim x|z\)。中介變數必須要是類別。

圖 4.14呈現每一年的時間趨勢。注意年度必須要是類別變數。我們可以調整以下參數：

• index.cond指定子圖的順序
• scales設定X軸或Y軸的刻度名稱跟位置
• strip指定每一個子圖裡面的背景顏色

xyplot(value ~ month | as.factor(year),
       data=tsp, type="l", col="black",
       index.cond=list(c(4,5,6,1,2,3)),
       scales = list(at=c(1,6, 12),
              labels=c("Jan.","Jun.","Dec.")),
       strip = strip.custom(bg="gray90"),
       ylab="percent")

Figure 4.14: 每年的意外死亡折線圖

本課程介紹的lattice套件繪圖功能，雖然品質不像ggplot2，但是提供ggplot2以外的選擇。

除了Deepayan Sakar的著作外，有興趣的同學可以參考Purdue大學的lattice教學資料，裡面很詳盡地說明各種圖形以及相關指令。

---

5 以ggplot 為基礎的套件

5.1 ggpubr

ggpubr 指的是ggplot publication ready。請下載ggpubr這個套件。詳細介紹可以參考這個網站。

ggpubr的優點是美化ggplot2的圖形，指令相當類似，例如圖5.1顯示鄉村與城市的死亡率長條圖，並且以虛線標示變數的平均值：

mor<-lattice::USRegionalMortality
ggpubr::gghistogram(mor, x = "Rate", bins=40,
   add = "mean", rug = TRUE,
   color = "Status", fill = "Status",
   palette = c("#EE0011", "#111100"))

Figure 5.1: ggpubr長條圖

圖5.2 則是箱型圖的應用，可以比較圖4.12，除了方向不同之外，還可以看到ggpubr的圖形顯示出觀察值。

p <- ggpubr::ggboxplot(mtcars, x = "cyl", y = "mpg",
          color = "cyl", palette =c("#00AFBB", "#E7B800", "#FC4E07"),
     add = "jitter", shape = "cyl")
 p

Figure 5.2: ggpubr的箱型圖

5.2 ggstatsplot

ggstatsplot套件結合圖形與統計，直接呈現圖形上面的統計結果。先下載套件：

r install.packages('ggstatsplot')

我們以ISLR::Default為例，進行交叉列表並卡方檢定，呈現圓餅圖 5.3：

# for reproducibility
set.seed(02139)
data(ISLR::Default)
# plot
ggstatsplot::ggpiestats(
  data = ISLR::Default,
  x = student,
  y = default,
  title = "Student and Default", # title for the plot
  legend.title = "student", # title for the legend
  caption = substitute(paste(italic("Source"), ": ISLR")),
  messages = FALSE
)

Figure 5.3: 卡方檢定圓餅圖

兩個連續變數之間可以用相關係數來描述兩者之間的關聯，而ggstatsplot可以顯示相關係數：

ggstatsplot::ggscatterstats(
  data = alr4::florida,
  x = Gore,
  y = Bush,
  xlab = "Vote for Gore",
  ylab = "Vote for Bush",
  title = "County-by-county vote for president in Florida in 2000 for Bush, Gore",
  caption = substitute(paste(italic("Source"), ": alr4")),
  messages = FALSE,
  bf.message = F,
  results.subtitle = T
)

Figure 5.4: ggstatsplot散佈圖與迴歸係數

6 函數與ggplot2

還記得一開始的課程目標是畫出類別變數與連續變數之間的關係，並且標示各類別的平均值？請見圖 6.1。

library(alr4)
data(salary)
p <- ggplot(salary, aes(x=rank, y=salary)) +
  geom_point() +
  theme_classic() 

p + ggplot2::stat_summary(fun='mean', colour='red',size=3, geom = "point")

Figure 6.1: ggplot2散佈圖加上平均值

我們運用stat_summary()計算各類別對應的變數的平均數，還可以加上標準差，並且運用pointrange的指令，顯示平均值與標準差連起來的直線，可以用來比較各類別之間的差異程度。圖 6.2顯示，教授的薪水明顯高於助理教授，但是不一定高於副教授。

p + stat_summary(fun = mean, colour='#2200ae',size=1.3,
               fun.min = function(x) mean(x) - sd(x), 
               fun.max = function(x) mean(x) + sd(x), 
               geom = "pointrange") +
  stat_summary(fun = mean, geom = "line") +
  scale_x_discrete(breaks=c('Asst','Assoc','Prof'),
    labels=c('Assistant Professor','Associate Professor', 'Full Professor'))

Figure 6.2: ggplot2散佈圖加上平均值及標準差

7 作業

7.1 基礎繪圖指令作業

1. 請針對mtcars這筆資料的wt以及mpg繪製散佈圖，並且加上車輛的名稱。

2. 請針對ChickWeight的Diet變數畫出直方圖，並且標示次數（提示：在圖形上面加上文字的語法為：text(x, y, labels)，x, y 分別是標示次數的位置，而標示次數要轉為字串）。

3. 請畫直方圖表示mtcars這筆資料中的變速箱（am）是屬於自動或是手動排檔之中，gear的比例。

4. 請嘗試讀取PP0797B2.sav這筆資料，然後畫長條圖表示partyid這個變數中，不同年齡(age)是否同意Q9（政府所做的是大多數是正確的）的相對次數。請注意讀取SPSS資料的方式，有的方式可以不讀取原先資料設定的遺漏值。

5. 請畫圖表示Nile的水位變化，並且指出大概是哪一年的水位最低與最高。

6. 請從CS3171D1A.csv資料選取前五個縣市，然後針對老年人口比率畫2000至2010年的折線圖，並且於圖例標示五個縣市。

7. 請畫箱型圖表示mtcars這筆資料中的gear的各個類別的hp分佈。

8. 請畫圖表示美國1986年調查的嬰兒體重(MASS::birthwt)的分佈。

9. 請用折線圖表示Orange資料裡五種樹的樹齡與樹圍之間的關係。（提示：請考慮用for迴圈加上五條折線以及五種線條與點狀。另外請用subset()指定每一種樹的資料。）

10. 請觀察reshape2::tips這筆資料的tips與total_bill兩個變數的關係，並且標記顧客是否抽煙(smoker)。請記得加上圖例。

7.2 ggplot2作業

1. 種族與社會經濟地位有關嗎？請用ggplot2繪製直方圖顯示hsb2的race與ses的關係。

2. 請針對flights這筆資料之中的arr_delay，畫圖顯示每一家航空公司(carrier)的分布情況：

3. 請嘗試篩選airquality這筆資料的8,9兩個月份的資料，並且用ggplot2畫出Ozone變數之折線圖。X軸是月份加上日期。如果有時間點沒有資料（遺漏值），請仍然保留該時間點。

4. 請參考Pew Research Center畫點狀圖顯示下列的資料：

library(dplyr)
DT <- tibble(Country=c("Spain", "France", "Germany","Italy", "Sweden"),
             Populist=c(52,58,65,50, 73),
             Mixed=c(70,72,70,65,81),
             Nonpopulist=c(85,84,88,68, 90))
kableExtra::kable_styling(knitr::kable(DT))

Country	Populist	Mixed	Nonpopulist
Spain	52	70	85
France	58	72	84
Germany	65	70	88
Italy	50	65	68
Sweden	73	81	90

提示：

• 請轉換為tidy data
  
• 請考慮使用ggthemes套件

5. 請畫出以下四筆資料的密度長條圖在同一張圖上（提示：運用reshape2套件中的melt以及lattice套件）

set.seed(2019)
m=1000; n=580; p1=0.54; p2=0.60; p3=0.70; p4=0.80
res1=rbinom(m, n, p1);res2=rbinom(m, n, p2); res3=rbinom(m,n,p3);
res4=rbinom(m, n, p4)

6. UCBAdmissions是一筆表格資料，表示每個系錄取與拒絕男性與女性的人數。請轉換這筆資料，然後用畫圖表示六個系男性與女性的錄取率，並回答哪一個系的男性與女性錄取率差異最大。（提示：使用dplyr指令以及dotplot()點狀圖）

7. 請畫圖顯示mtcars這筆資料中，汽車是否為automatic與mpg的關係，並且用文字說明兩種變速箱的汽車的mpg集中趨勢。

8. 請用UsingR的baycheck資料，畫圖表示1960到1986年的Bay Checkerspot蝴蝶的數量變化：

9. 請根據ISLR::College這筆資料的Top10perc排序，篩選出新生為高中前10%畢業生的比率最高的前15名學校，觀察錄取率(Accept/Apps)與生師比(S.F. Ratio)的關係，並且顯示學校的名稱。

10. 請把Tondutrend.csv這筆資料轉換成時間序列資料，並且用xyplot分別畫四個態度的趨勢。請注意X軸的刻度以及子圖的順序。

8 更新日期

最後更新日期 03/26/2021