1 課程目標
2 基礎繪圖指令
3 ggplot2繪圖
4 qplot繪圖
5 lattice 繪圖
6 以ggplot 為基礎的套件
- 6.1 ggpubr
- 6.2 ggstatsplot
7 函數與ggplot2
8 作業
- 8.1 基礎繪圖指令作業
- 8.2 ggplot2以及lattice作業
9 更新講義時間

1 課程目標

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

file<-here("data","studentsfull.txt")
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")

Figure 1.1: 類別平均值的散佈圖

又例如用Emoji做折線圖1.2：

tmp1 <- nycflights13::weather
tmp1 <- tmp1[which(tmp1$origin=='JFK'), ]
 tmp2 <- tmp1 %>% 
     group_by(month) %>%
      summarise(Temp=mean(temp))


p1 <- ggplot(aes(x=as.factor(month), y=Temp, group=1), data = tmp2) +
   geom_line(stat = 'identity') +
   labs(x = "Month", y = "Temperature")

p1  + geom_emoji(emoji = "1f31e")

Figure 1.2: Emoji的折線圖

在使用表情符號之前必須先下載emoGG套件，語法如下：

# downloading the package from GitHub
devtools::install_github(
  repo = "dill/emoGG", # package path on GitHub
  dependencies = TRUE,                # assumes you have already installed needed packages
  quick = TRUE                         # skips docs, demos, and vignettes)

我們還可以透過資料視覺化展現不同團體對於同一政策的立場如圖1.3：

#read data
file<-here::here('data','PP2197E5C.sav')
dat<-sjlabelled::read_spss(file)

library(reshape2)
dat$bluegreen<-car::recode(dat$partyid, "1=1;2=2;3=1;
                           4=1;5=2;6=2;7=3;8:9=4")
#function
myrecode <- function(x){
recode(x, '1'=0, '2'=0, '3'=1, '4' = 1)
}

dat<-dat %>% mutate_at(c("Q1","Q2","Q3","Q5","Q7","Q10", "Q11","Q15","Q18","Q19","Q20","Q21","Q22","Q23",
                    "Q24","Q25"), myrecode)

varname<-c(1,2,3,5,7,10,11,15,18,19,20,21,22,23,24,25)
S<-paste('Q',varname, sep='')

listtable<-lapply(dat[,S], function(x)   xtabs(~x + dat$bluegreen))

MyList <- lapply(listtable, prop.table, 2)

k <- matrix(rep(NA,64),nrow=16,ncol=4)
for (i in 1:16){
  k[i,]<-as.vector(MyList[[i]][2,])
}
k<-as.data.frame(k)
#data set
rownames(k)<-NULL
colnames(k)<-c("Pan-blue","Pan-green","TPP","IND")
#reverse the order of questions so Q1 changes to Q25
k$Q<-noquote(paste('Q', rev(varname), sep=''))

tmpo<-melt(k, id.vars=c("Q"))

tmpo$Qnew <-as.character(tmpo$Q)
#tmpop$value<-100*tmpop$value
tmpo<-tmpo%>% mutate(party=variable,
                       value=100*value,
                       Qnew=factor(Qnew,levels=S))
showtext::showtext_auto()
ylabel<-c('滿意蔡英文','滿意蘇貞昌','滿意陳時中','滿意紓困措施','校正回歸適合','對指揮中心信心','對政府防疫信心','地方可以買疫苗','願打中國疫苗','願打代理歐美疫苗','滿意買疫苗努力','指揮中心決定正確',
               '指揮中心考慮民眾福利',
        '相信陳時中','政府可以用個資','擔心資安被侵害')
G1<-ggplot(data=tmpo, aes(x=Qnew,y=value,label=value)) +
  geom_point(stat="identity", aes(color=party), size = 6) +geom_text(aes(label=sprintf("%0.1f", value)),vjust = 2, size=3) +
  labs(caption='資料來源：請見中華民國政大選研之友協會臉書粉絲頁version 4',
  title='', x='', y='') +
   geom_vline(xintercept = 50, linetype='dashed') +
   theme_economist() +
theme(text=element_text(size=14,family='Georgia'),
      plot.caption = element_text(size=11),
  axis.text.x = element_text(size=14, family='YouYuan'),
  plot.title = element_text(family='YouYuan', size=10),
        legend.title = element_blank())   +
  guides(fill=FALSE) +
  coord_flip()+
  scale_color_manual(values=c('#3322EE','#2AA111','white','yellow','#CDC1B2', '#B1000B'))+
  scale_x_discrete(breaks=S,
   labels=rev(ylabel))
G1

Figure 1.3: 民眾對疫情看法

2 基礎繪圖指令

有人說，「一張照片勝過千言萬語」。越來越來媒體報導、研究報告、論文等等，以統計圖形表現變數之間的關係，讓讀者可以很快地了解文字內容，這就是所謂「資料視覺化」。當然，統計圖形也可能有誤導作用，必須要小心解讀。

R內建了許多基本的繪圖指令，讓我們可以很輕鬆地畫出專業的圖形，但是ggplot2也有許多技巧，所以我們必須學習並且熟練這些語法。一旦熟悉之後，你可以客製化你想要的感覺。

2.1 基本元素

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

col: 控制符號的顏色，也可以是灰階的濃度，注意不要打成color或者colour。
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: 控制圖形中符號的大小
cex.axis: 控制X或是Y軸的數字刻度大小
cex.main: 控制圖形標題的大小
cex.lab: 控制X或是Y軸標題的文字的大小
cex.names: 控制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軸任意的刻度。

2.2 直方圖(Bar Plot)

直方圖可以表現單一類別變數的分佈，讓人一眼就可以看出哪一個類別有最多的次數或者百分比，但是不管是次數或者比例，必須要先把原始資料轉成表格。例如我們想呈現學生的科系的次數分佈，必須先轉成表格如表 2.1：

TXT<-here("data","studentsfull.txt")
students<-read.table(TXT, header=T, sep='')
stu <- table(students$Department)
kableExtra::kable_styling(knitr::kable(stu, caption="科系分佈的次數"))

Table 2.1: 科系分佈的次數
Var1	Freq
Aerospace	4
Chemistry	3
Economics	4
English	3
Journalism	4
Mechanics	5
Physics	3

直方圖的對象應該是具有名稱的摘要，所以要先把數字、字串或者是類別的向量改為表格的形式。再用barplot()產生直方圖 2.3，顯示學生的科系的次數分佈：

barplot(stu, main="Departments",  xlab="", ylab="frequency", 
        cex.names = 0.4)

Figure 2.3: 直方圖一

再舉一個例子，我們用diamonds這筆資料的cut變數，顯示其直方圖 2.4：

o.par <- par()
par(bg = '#a0a1c2')
barplot(table(diamonds$cut))

Figure 2.4: 鑽石切割等級的直方圖

par(o.par)

直方圖也可以指定直方本身的顏色以及外框的顏色，例如圖 2.5 換了新的顏色。另外注意我們用cex.names控制X軸類別的文字大小。

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

Figure 2.5: 直方圖二

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

如果圖例所使用的變數是類別變數，圖例函數寫成： 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'), cex.axis= 0.5,
        legend = levels(students$Gender))

Figure 2.6: 直方圖三

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

圖例的位置有“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'), cex.axis=0.7)
legend("top", fill=c('blue', '#EE330011'), c("M","F"))

Figure 2.7: 直方圖與圖例

我們也可以設定圖的外框：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.8: 圖例在圖外面

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

Figure 2.9: 圖例在圖外面

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

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

Figure 2.10: 直方圖四

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

2.3 盒型圖 (Box Plot)

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

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

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

Figure 2.11: 箱型圖一

中間的粗線代表中位數。盒型圖上方的線稱為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.12 對照上述的計算結果：

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.12: 盒型圖二

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

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

Figure 2.13: 盒型圖三

然後請輸入下列指令：

y<-state.abb

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

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

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

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

Figure 2.14: 盒型圖四

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

2.3.1 課間練習

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

2.4 長條圖 (Histogram)

我們可以從長條圖觀察連續變數的分佈，例如偏態以及峰值等等。因為連續變數有一定數量的值，而非像類別變數的值是離散的，所以應該用長條圖來呈現。例如圖 2.15 顯示美國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.15: 長條圖一

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

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

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

Figure 2.16: 長條圖二

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

我們用模擬的資料來表現長條圖的參數break的用途。我們從常態分佈的x抽出100個，計算平均值後記錄為y，然後重複1000次，也就是有1000個x的平均值。在圖 2.17，左邊的圖最多有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.17: 兩個長條圖並列

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

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.18: 加上機率密度曲線的長條圖

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

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

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

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

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.19: 鑽石價格機率密度

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

2.5 折線圖 (Line Chart)

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

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

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

Figure 2.20: 空氣品質折線圖一

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

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

Figure 2.21: 空氣品質折線圖二

2.5.1 時間序列資料

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

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

Figure 2.22: Lake Huron折線圖1

可以下載tsbox這個套件，針對ts的物件，畫出折線圖：

library(tsbox)
ts_plot(LakeHuron)

Figure 2.23: Lake Huron折線圖2

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

回到airquality這筆資料。airquality 是時間序列資料，應該要用plot.ts(ts())來畫折線圖。我們轉換Wind這個變數為時間序列，注意frequency = 31

x <- as.vector(airquality$Wind)
my_ts <- ts(x, frequency = 31, start=c(5,1),
            end=c(9,30))
plot.ts(my_ts, type='b')

Figure 2.24: 風速趨勢圖

圖 2.24顯示5, 6, 7, 8, 9這五個月的風力。
第二個例子是用我們自己收集的統獨趨勢資料，裡面有統一、獨立、維持現狀、無反應等四個變數，也就是要畫四條線。首先轉換資料：

CSV<-here::here("data","Tondutrend.csv")
trend<-read.csv(CSV, header=T, sep=",")
tonduts<-ts(trend, start=c(1992.1), frequency=2)
tonduts
## Time Series:
## Start = 1992.1 
## End = 2020.1 
## Frequency = 2 
##       統一 維持現狀  獨立  拒答
## 1992  2.40    66.10  7.70 23.80
## 1993 18.60    42.20  4.00 35.20
## 1993 20.40    51.40 11.40 16.70
## 1994    NA       NA    NA    NA
## 1994    NA       NA    NA    NA
## 1995    NA       NA    NA    NA
## 1995 19.10    46.10 11.10 23.80
## 1996 22.30    41.10 12.20 24.40
## 1996 21.00    48.50 14.10 16.40
## 1997 27.80    47.10 16.10  9.00
## 1997 22.80    50.80 17.90  9.00
## 1998 20.50    49.90 17.40 12.20
## 1998 20.30    50.50 20.00  9.20
## 1999 18.80    53.80 18.90  8.50
## 1999 13.30    58.80 17.80 10.10
## 2000 20.80    52.50 16.20 10.60
## 2000 26.80    51.10  9.00 13.00
## 2001 23.80    55.50 13.30  7.40
## 2001 26.20    52.50 15.40  6.00
## 2002 21.50    53.50 17.90  7.10
## 2002 21.60    52.90 18.80  6.80
## 2003 17.60    56.60 19.20  6.60
## 2003 14.90    55.90 21.60  7.70
## 2004 12.80    60.20 19.60  7.40
## 2004 14.10    57.70 20.70  7.50
## 2005 15.70    59.30 19.80  5.20
## 2005 16.20    58.40 20.10  5.20
## 2006 14.90    60.90 18.50  5.70
## 2006 15.90    60.70 18.60  4.80
## 2007 18.20    53.30 23.10  5.40
## 2007 12.90    58.20 20.80  8.00
## 2008 11.50    59.60 20.40  8.50
## 2008 10.50    59.60 25.70  4.10
## 2009 10.20    63.10 21.30  5.50
## 2009 10.60    63.90 19.60  5.90
## 2010 12.30    61.00 22.80  3.80
## 2010 10.00    64.60 21.00  4.40
## 2011 11.10    61.30 23.00  4.60
## 2011 10.80    63.70 19.10  6.40
## 2012 10.70    64.80 19.10  5.40
## 2012 11.10    62.80 19.60  6.50
## 2013 12.30    59.20 23.30  5.20
## 2013 11.60    61.30 21.70  5.40
## 2014 10.90    60.80 22.90  5.40
## 2014  8.80    61.60 23.60  5.90
## 2015  9.60    61.60 20.50  8.40
## 2015 10.80    60.60 22.20  6.40
## 2016 10.00    61.30 22.80  5.90
## 2016 11.40    60.10 22.50  5.90
## 2017 13.20    58.70 23.10  5.00
## 2017 14.40    59.80 19.90  5.90
## 2018 16.50    58.10 20.00  5.50
## 2018 17.50    57.90 19.50  5.20
## 2019 10.60    57.60 25.80  6.00
## 2019  7.50    57.60 28.00  6.98
## 2020  8.61    56.56 30.30  4.54
## 2020  7.83    56.08 30.96  5.13

相對於其他軟體，R的優點在於允許時間點有缺漏，畫出來的圖形就有缺少某個或是某些時間點，不會強制連接前後兩個時間點。接下來是畫折線圖。圖 2.23 顯示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(2,3,1,2),ylab="%",xlab=NULL,pch='1', lwd=3,frame.plot=F,col=c("gray20","gray60", "black", "gray80"),xaxt="n", main="台灣民眾的統獨立場, 1992.6-2021.12")
axis(1, at=seq(1992,2021,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(2,3), bty='n', lwd=3)
legend("bottomright", c("獨立","無反應"), inset=c(0, -0.6), col=c("black", "gray80"),lty=c(1,2), bty='n', lwd=3)
text(2000, 30, paste("第一次政黨輪替"))

Figure 2.25: 統獨趨勢折線圖1

在圖形 2.25中，我們設定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這個指令。相信這個折線圖會讓你的報告比其他軟體畫的圖看起來更專業。

我們也可以用tsbox套件中的ts_plot()來畫圖??：

par(family='TaipeiSansTCBeta')
png('./Fig/unplot.png')
ts_plot(tonduts, title='')
dev.off()

gra1<-here::here('Fig','unplot.png')
knitr::include_graphics(gra1, dpi=100)

2.5.2 顯示中文字

這張圖裡面有中文，所以練習時，我們需要設定中文字型，不然會顯示亂碼。如果你打開字體簿，可以找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.3 課間練習

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

2.6 散佈圖 (Scatter Plot)

接續一開始的圖例說明，我們遇到兩個變數時，通常用散佈圖來顯示兩者之間的關係，指令就是plot()。例如圖 2.26：

with(anscombe, plot(x1, y1, pch=16, cex=1.2, col='gray25'))

Figure 2.26: 散佈圖

我們可以先幫觀察值畫一個框，再用points這個函數，如圖2.27：

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.27: 散佈圖二

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

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.28: 散佈圖三

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

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

也可以用text貼上對應觀察值的數字如圖 2.30：

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.30: 散佈圖五

我們用state.x77這筆資料再來示範一次帶有資料名稱的散佈圖。首先把這個矩陣改為資料框，然後選兩個變數畫圖，最後加上各州的州名，州名來自於rownames()這筆資料的每一列的名稱。

x77 <- data.frame(state.x77)
with(x77, plot(Income, HS.Grad))
with(x77, text(Income, HS.Grad, rownames(x77), adj = -0.2, col='blue'))

Figure 2.31: 各州收入與高中畢業比例散佈圖

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

2.7 特殊點狀圖

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

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

Figure 2.32: 特殊點狀圖一

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

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

Figure 2.33: 特殊點狀圖二

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

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.7.1 課間練習

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

2.8 輸出圖形

我們可以輸出圖形到指定的路徑，方便放在報告當中。例如我們產生一個圓餅圖，然後輸出為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, fill,…))+ geom_object()+theme()

aes()裡面指定X與Y軸的變數，以及用group指定第三個變數，用fill指定塗色的變數。
geom_object可以是geom_bar（直方圖）、geom_point（散佈圖）、geom_line（折線圖）、geom_histogram（長條圖）等等。
geom_object可設定圖形的顏色、形狀、大小等等，例如geom_bar(fill=‘red’)，也可以設定顏色、形狀、大小隨著特定變數改變。
繪圖的對象可以是資料框或是已經轉為表格的現成資料。
繪圖之前應該特別注意變數的型態。
有興趣的同學可參考這個網頁，上面有許多ggplot2的例子。

3.1 直方圖 (Bar plot)

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

ggplot(ChickWeight, aes(x=Diet)) +
  geom_bar(fill='#ee22ffaa') +
  ggtitle('Diet')

Figure 3.1: ggplot2直方圖1

我們可以在ggplot告訴系統x指的是什麼變數，然後使用 bar 這個 geom呈現圖形3.1。在這個例子，我們要呈現的是一個類別變數的次數分配。
在產生直方圖時，我們也可以使用 stat=‘count’ 這個統計轉換來計算直方圖每個 bin 的數值，然後再使用 bar 這個 geom 呈現類別變數的圖形3.2。

ggplot(ChickWeight, aes(x=Diet)) +
  geom_bar(stat='count', fill='#11ff44cc') +
  ggtitle('Diet')

Figure 3.2: ggplot2直方圖2

或者是先轉換成表格，但是要設stat=“identity”，畫成3.3：

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="#44bb33") + ylab('Count')
g1

Figure 3.3: ggplot2直方圖3

stat=‘identity’用在對x與y繪圖的時間，也就是科系是x，人數是y，直方圖呈現每個科系的人數。
也可以直接使用原始資料students畫圖3.4：

ggplot(students, aes(x=Department)) + 
  geom_bar(fill='#33aacc')

Figure 3.4: ggplot2直方圖4

我們也可以自己輸入資料，用現成的次數分配畫圖。先輸入四家公司的市值：

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

然後我們轉成資料框，再畫圖 3.5。我們可以用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.5: 次數分配轉成直方圖

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

3.1.1 課間練習

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

3.1.2 加上百分比到Y軸

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

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.6: 加上百分比的直方圖

3.1.3 加上數據到直方圖

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

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

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

3.1.4 轉成水平

我們可以轉成水平，加上百分比的標籤：

# percentage
stu.t <- stu.t %>% mutate(prop.d=round(Freq/sum(Freq),3))

g2 <- ggplot(stu.t, aes(x = Dep, y=prop.d, fill=Dep)) +
  geom_bar(stat="identity", position="dodge",
            color="black", width = 0.75)+
 geom_text(aes(x=Dep, y=prop.d,
  label=sprintf("%1.1f%%", 100*prop.d)), size=8, colour='gray40', hjust=1) +
  scale_y_continuous(labels = scales::percent)  +
  labs(y='%') +
  coord_flip()  +
  theme(legend.position = 'none')

g2

Figure 3.8: 加上資料標籤的水平直方圖

3.1.5 應用stat_count()

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

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.9: ggplot2直方圖三

用geom_bar()得到同樣的直方圖：

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

Figure 3.10: ggplot2直方圖四

我們接下來用stat(prop)計算每一個類別所佔的全部比例，並且表現在y軸。然後用stat_count加入每一個類別的比例數字如圖3.11：

#ggplot2
p1<-ggplot(salary)+geom_bar(aes(x=rank, y=stat(prop), group=1), fill='#bb00cc')
#layer
p1+stat_count(aes(x=rank,y=stat(prop), group =1, label=round(100*stat(prop),digits = 3)), 
      geom='text', color='white',vjust=1.3) +
  scale_y_continuous(label = scales::percent)

Figure 3.11: ggplot2直方圖五

3.1.6 兩個類別變數的直方圖

我們想用直方圖表現兩個變數之間的關係，例如在學生成績這個資料中，我們想知道每一個系考試超過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,caption ='學生成績資料'))

Table 3.1: 學生成績資料
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.12：

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

Figure 3.12: 並列直方圖

如果要讓色塊表示條件機率，我們可以用 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, caption='學生成績統計'))

Table 3.2: 學生成績統計
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

比較表 3.2 跟表3.1，可以發現前者是後者的統計，也就是統計有幾個人通過，例如航太系有4人通過，化學系有2人通過，1人沒通過等等。

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

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

Figure 3.13: ggbarplot2堆疊直方圖

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

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

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.14: ggplot2堆疊直方圖

3.2 用tabyl製作直方圖

3.2.1 單一變數

janitor套件中有一個產生次數分配的函式：tabyl，透過%>%這個符號可以結合ggplot2，例如：

t1.g <- students %>%
  janitor::tabyl(Department) %>%
  ggplot() + aes(x=Department,y=n, group=1) + 
  geom_bar(stat='identity', fill='#ABC013') 
t1.g

Figure 3.15: tabyl直方圖ㄧ

在以上的語法中，我們先產生一個表格，裡面有各個科系的人數以及百分比，如表3.3。我們就用這個表格當作ggplot2的資料，產出圖3.15。

t1 <- students %>%
  janitor::tabyl(Department) 
kable(t1, caption='學生系別次數分配')

Table 3.3: 學生系別次數分配
Department	n	percent
Aerospace	4	0.1538
Chemistry	3	0.1154
Economics	4	0.1538
English	3	0.1154
Journalism	4	0.1538
Mechanics	5	0.1923
Physics	3	0.1154

圖3.16則是應用百分比，但是要把Y軸轉變為百分比。

t2 <- students %>%
  janitor::tabyl(Department) %>%
  ggplot() + aes(x=Department,y=percent, group=1) + 
  geom_bar(stat='identity', fill='#11FF00AA') +
  scale_y_continuous(label = scales::percent) 
t2

Figure 3.16: tabyl直方圖二

3.2.2 兩個變數

接下來我們嘗試畫兩個變數交叉表的圖形。在開始前，先用tabyl產出交叉表 3.4：

t3 <- students %>%
janitor::tabyl(Department, Gender) %>% 
  janitor::adorn_percentages('col') 
kableExtra::kable_styling(knitr::kable(t3, caption='性別與科系交叉表'))

Table 3.4: 性別與科系交叉表
Department	F	M
Aerospace	0.2308	0.0769
Chemistry	0.0769	0.1538
Economics	0.2308	0.0769
English	0.0000	0.2308
Journalism	0.2308	0.0769
Mechanics	0.1538	0.2308
Physics	0.0769	0.1538

觀察這個表格，發現資料在兩個欄位，所以我們用melt()把資料放在一個欄位，也就是科系會重複出現兩次，先排女生(F)的科系比例再排男生(M)的科系比例：

t3.n <- reshape2::melt(t3)
kable(t3.n)

Department	variable	value
Aerospace	F	0.2308
Chemistry	F	0.0769
Economics	F	0.2308
English	F	0.0000
Journalism	F	0.2308
Mechanics	F	0.1538
Physics	F	0.0769
Aerospace	M	0.0769
Chemistry	M	0.1538
Economics	M	0.0769
English	M	0.2308
Journalism	M	0.0769
Mechanics	M	0.2308
Physics	M	0.1538

把以上的步驟結合ggplot2得到以下的堆疊圖3.17：

t3.g <- students %>%
janitor::tabyl(Department, Gender) %>% 
  janitor::adorn_percentages('col')  %>%
  melt() %>%
  ggplot() + aes(x=Department,y=value, fill=variable) + 
  geom_bar(stat='identity') +
  scale_y_continuous(label = scales::percent) 
t3.g

Figure 3.17: tabyl直方圖三

比較圖3.13以及3.17，如果我們要知道每一個系的性別比例，還要加工一下。我們先產出帶有次數以及百分比的交叉表，確定百分比正確：

t4 <- students %>%
     janitor::tabyl(Department, Gender) %>% 
     janitor::adorn_percentages('row') %>% adorn_ns()
kableExtra::kable_styling(knitr::kable(t4, caption='性別與系別交叉表（系為100%）'))

Table 3.5: 性別與系別交叉表（系為100%）
Department	F	M
Aerospace	0.7500 (3)	0.2500 (1)
Chemistry	0.3333 (1)	0.6667 (2)
Economics	0.7500 (3)	0.2500 (1)
English	0.0000 (0)	1.0000 (3)
Journalism	0.7500 (3)	0.2500 (1)
Mechanics	0.4000 (2)	0.6000 (3)
Physics	0.3333 (1)	0.6667 (2)

我們把所有步驟透過%>%串連，送到ggplot2，就大功告成，見圖3.18:

t4.g <- students %>%
     janitor::tabyl(Department, Gender) %>% 
     janitor::adorn_percentages('row') %>%
  melt() %>% 
  ggplot()+
  aes(x=Department,y=value,fill=variable) + 
  geom_bar(stat='identity') +
  scale_y_continuous(label = scales::percent) +
  scale_fill_manual(values=c('#C4106A','#12026C'))
t4.g

Figure 3.18: tabyl直方圖四

3.2.3 課間練習

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

3.3 長條圖 (Histogram)

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

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

Figure 3.19: ggplot2長條圖一

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

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

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

Figure 3.20: ggplot2長條圖二

3.3.1 課間練習

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

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

3.4 散佈圖 (scatter plot)

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

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

Figure 3.21: ggplot2散佈圖一

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

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.22: ggplot2散佈圖二

3.4.1 加上迴歸線

散佈圖可以用geom_smooth加上迴歸線表現兩個變數的線性相關，我們用method=‘lm’設定。如圖3.23：

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

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

也可以用非線性迴歸線表示兩者的相關，，我們用method=‘loess’設定。loess 代表Local Polynomial Regression. 如圖3.24：

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

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

3.4.2 課間練習

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

3.5 三個變數之散佈圖

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

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

Figure 3.25: 三個變數散佈圖一

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

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.26: 三個變數散佈圖二

也可以隨著另一個連續變數而改變顏色深淺如圖 3.27：

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

Figure 3.27: 三個變數散佈圖三

3.5.1 課間練習

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

3.6 盒型圖 (Box plot)

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

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

Figure 3.28: 盒型圖一

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

g1 <-ggplot(mtcars, aes(x=as.factor(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.29: 盒型圖二

由於 cyl 是數字變數，所以圖 3.28以及圖 3.29用as.factor把數字轉為因素。
當我們需要改變在X軸的類別變數的名稱時，可以用 scale_x_discrete(breaks=…,labels=c(…))加以調整，如圖 3.30：

o1<-ggplot(Orange, aes(x=Tree, y=circumference)) +
  geom_boxplot(aes(colour=Tree), size=1.2) 
o1+scale_x_discrete(breaks=1:5,
     labels=c("A","B", "C", "D", "E"))

Figure 3.30: 橘子樹箱型圖1

對照圖3.30的圖例與X軸，可以發現R會根據原始的類別數字1,…5，改為A,…E。我們如果希望按照數字順序排列，就先增加一個類別變數，順序改為1到5，然後再改為A到E，如圖3.31：

Orange <- Orange %>% mutate(trees=factor(Tree,
          levels=c(1,2,3,4,5)))
o1<-ggplot(Orange, aes(x=trees, y=circumference)) +
  geom_boxplot(aes(colour=trees), size=1.2) 
o1+scale_x_discrete(breaks=1:5,
     labels=c("A","B", "C", "D", "E"))

Figure 3.31: 橘子樹箱型圖2

如果我們只想看到數值落在特定範圍的類別，例如我們只想顯示50到150的馬力，超過200以上的資料就會被去掉，可以用limits加以控制，例如圖 3.32：

g1 + scale_y_continuous(limits = c(50, 150)) +
              theme_bw()

Figure 3.32: 盒型圖四

3.6.1 課間練習

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

3.7 甜甜圈圖

基礎指令可以畫圓餅圖，甜甜圈圖同樣可以顯示資料的次數分佈。例如我們想知道ggplot2裡面的的diamonds資料中的cut的次數分配，先針對這個變數計算個數如表 3.6：

dt <-group_by(diamonds, cut)
dat = summarize(dt, count=n())
kableExtra::kable_styling(knitr::kable(dat, caption='鑽石切割等級次數分佈'))

Table 3.6: 鑽石切割等級次數分佈
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.33：

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="", x="", y='')
p1

Figure 3.33: 甜甜圈圖一

xmax不能太大，不然畫圖會失敗。xmin越小，環裡面的面積越小，當xmin=0時，會得到圓餅圖3.34：

p1 = ggplot(dat, aes(fill=cut, ymax=ymax, ymin=ymin, xmax=4, xmin=0)) +
     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="Diamond Cut",x="",y='')
p1

Figure 3.34: 甜甜圈圖二

3.7.1 課間練習

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

3.8 折線圖

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

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

Figure 3.35: ggplot折線圖一

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

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

Figure 3.36: ggplot折線圖二

應用paste0()把原來的月份與日變數轉換成日期：

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

然後用日期畫折線圖 3.37：

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

Figure 3.37: 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 2024-05-01     Wind   7.4
## 2 2024-05-02     Wind   8.0
## 3 2024-05-03     Wind  12.6
## 4 2024-05-04     Wind  11.5
## 5 2024-05-05     Wind  14.3
## 6 2024-05-06     Wind  14.9

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

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

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

Figure 3.38: ggplot2折線圖四

3.8.1 課間練習

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

3.9 特殊圓點圖

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.39：

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.39: ggplot圓點圖一

然後加上代表係數的圓點，並加上係數的名稱，以及把X軸與Y軸對調，就大功告成。圖 3.40 顯示是否自動排檔對於依變數的影響相對來得大，但是有可能等於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.40: ggplot圓點圖二

3.9.1 課間練習

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

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.41: ggplot2係數圖

3.10 特殊風格

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

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.42: 主題風繪圖

3.11 結合統計與圖形

之前我們的圖形顯示資料的原始測量結果，但是我們可以在圖裡面直接呈現變數的統計值，例如盒型圖 3.29。這裡我們用直方圖 3.43 比較不同切割等級的鑽石的平均價格：

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

Figure 3.43: 平均值直方圖

我們以pscl這個套件裡面的admit資料為例，用直方圖 3.44顯示不同分數等級的人數：

d1 <- ggplot(pscl::admit, aes(score)) 
d1 + geom_bar(fill = '#b4b106')

Figure 3.44: 分數等級之直方圖

我們再用直方圖 3.45顯示不同分數等級的平均字彙分數。注意到這裡用stat_summary來顯示統計值：

d1 + stat_summary(aes(y = gre.verbal), fun = 'mean', geom='bar', fill='#3bb010')

Figure 3.45: 分數等級之平均值直方圖

我們以ggplot2這個套件裡面的mpg資料為例，圖 3.46顯示不同車型在市區行駛的每加侖哩程數，這裡我們用stat_summary來顯示統計值：

ggplot(mpg, aes(class, cty)) +
  geom_point() +
  stat_summary(geom = 'pointrange',  
               fun.y = mean, colour = "red", size = 0.6)

Figure 3.46: 平均值點狀圖ㄧ

圖 3.47則是同時顯示中位數以及平均值，方法是在fun分別填入median與mean。

ggplot(mpg, aes(class, cty)) +
  geom_point() +
  stat_summary(fun='mean',  colour = "orange", size = 0.6) +
  stat_summary(fun='median',  colour = "blue", size = 0.6)

Figure 3.47: 平均值點狀圖二

我們也可以在箱型圖加入一條橫線表示平均值，方式是加入一個對應x產生y平均值的函數（fun.data），例如我們想知道三種驅動方式（前輪、後輪、四輪）的市區每加侖的平均哩程數以及中位數，可畫圖如3.48：

G1<-ggplot(mpg, aes(as.factor(drv), cty, fill = as.factor(drv))) +
  geom_boxplot() +
  stat_summary(geom = 'crossbar',  
    width=0.8, fatten=0.5, colour='white',
      fun.data=function(x){return(c(y=mean(x),
            ymin=mean(x),ymax=mean(x)))}) +
             theme(legend.position = 'none')
G1

Figure 3.48: 平均值箱型圖

3.12 Likert圖

Likert是一種階層式量表，例如政治容忍量表，通常包括三題，例如容不容忍主張激進主義的學者出版相關書籍、到校園演講、到電視上發表言論等等，第一題的容忍度最低，第三題的容忍度最高，容忍度越低，越不會同意後面的題目。
HH套件裡面有一個函式叫做$\texttt{likert}$，可以用來分析多個題目的同意程度，呈現每一題同意的多寡，讓讀者比較哪些題目有相對比較高的同意程度。
應用這個函式分析不同背景民眾在一個問題上的同意程度，例如性別、教育程度、年齡等類別在一個立場上的差異，如圖3.49：

library(HH)
file <- here::here('data','HXC23014 Harvard Poll Data.sav')
dt <- sjlabelled::read_spss(file)
dt$pid<-c()
dt$pid[dt$Q61==1] <- 1
dt$pid[dt$Q61==2] <- 2
dt$pid[dt$Q61==3 & dt$Q62==1] <- 2
dt$pid[dt$Q61==4 & dt$Q62==1] <- 2
dt$pid[dt$Q61==5 & dt$Q62==1] <- 2
dt$pid[dt$Q61==3 & dt$Q62==2] <- 1
dt$pid[dt$Q61==4 & dt$Q62==2] <- 1
dt$pid[dt$Q61==5 & dt$Q62==2] <- 1
dt <- dt|>mutate(Gender = recode_factor(Q1, `1` = 'Male', `2`='Female', 
                                        .default = NA_character_)) |>
  mutate(heritage = dplyr::recode_factor(Q32,
                                        `1`="Heritage Doesn't Influence",
                                        `2`="Heritage Have Influence")) |>
  mutate(pid = dplyr::recode(pid, 
                             `1`='Democrats', `2`='Republicans')) |>
  mutate(nD103 = as.numeric(D103)) |>
  mutate(nD103 = car::recode(nD103, "1:2=1;3=2;4:5=3;6:8=4")) |>  
  mutate(edu = dplyr::recode_factor(nD103, 
  `1`="Less than high school",
  `2`="High school",
  `3` = "Some college",
  `4` = "Bachlor degree")) |>
  mutate(age = car::recode(Q57, "18:34=1;35:49=2;50:64=3;65:90=4;95:98=NA")) |>
  mutate(age.f = recode_factor(age, `1`= "18-34 years old",
                               `2`= "35-49 years old",
                               `3`= "50-64 years old",
                               `4`= "65-90 years old",
                               .default = NA_character_))

likertdf24 <- dt %>% dplyr::select(heritage, edu, age.f,  Gender) 

t0<-with(likertdf24, table(heritage))
t1<-with(likertdf24, table(Gender, heritage))
t2<-with(likertdf24, table(age.f, heritage))
t3<-with(likertdf24, table(edu, heritage))

t.all <- data.frame(rbind(t0, t1, t2, t3))
t.all

                  Heritage.Doesn.t.Influence Heritage.Have.Influence

t0 924 723 Male 414 401 Female 506 318 18-34 years old 231 279 35-49 years old 273 282 50-64 years old 242 97 65-90 years old 177 64 Less than high school 59 44 High school 257 113 Some college 301 150 Bachlor degree 307 416

colnames(t.all)<-c("Heritage Doesn't Influence", "Heritage Have Influence")

alltmp1 <- data.frame(Question=row.names(t.all), t.all[,1:2])
rownames(alltmp1) <- c(1:11)
alltmp1$Subtable <- c('Heritage Have Something with You?',
                      rep('Gender',2), rep('Age', 4),
                      rep('Education', 4))

colnames(alltmp1)<-c('Question',"Heritage Doesn't Influence", "Heritage Have Influence",'Subtable')
alltmp1[1,1]<-'All survey responses'

p1 <- HH::likert(Question ~.|Subtable, alltmp1, as.percent = TRUE, 
           positive.order = F, 
           scales = list(y = list(relation = "free")), layout = c(1, 5),
           main='', ylab='')

p1

Figure 3.49: 傳統態度與人口背景

4 qplot繪圖

qplot是ggplot2的簡易版，但是也可以產生高品質的圖形，例如我們想畫兩個隨機亂數的散佈圖（4.1）：

x<-rbinom(1000, 610, 0.2)
y<-rgamma(1000, 3)
qplot(x, y, color='#4b4102') +
  theme(legend.position = 'none')

Figure 4.1: qplot散佈圖1

可以看出不需要變數在資料框，就可以畫圖。而且語法更為直觀，但是可以加入ggplot2的控制設定。
利用geom加入折線，例如圖4.2：

x <- seq(-2, 2, length = 100)
y <- x^2
qplot(x, y, geom = "line")

Figure 4.2: qplot函數圖

4.1 qplot 直方圖

我們接下來畫單一變數的直方圖，用geom = c(‘bar’)：

library(carData)
qplot(data=WeightLoss, group, y = ..count../sum(..count..), 
      fill = group, geom = c('bar')) +
   scale_fill_brewer(palette = "Set5") +
  scale_y_continuous(labels = scales::percent) +
  labs(y='')

Figure 4.3: qplot直方圖

4.2 qplot 直方堆疊圖

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

library(ggplot2)
qplot(cut, data = diamonds, 
      geom = c("bar"),   fill = color, group = color)

Figure 4.4: qplot直方堆疊圖

如果用ggplot，先以N=n()計算每個成色在每個類別中的原始數目，然後把N放在Y軸，切割等級放在X軸。請見圖 4.5。

library(ggplot2)
df <- diamonds %>% dplyr::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 4.5: 灰階原始數目直方圖

4.3 qplot散佈圖

qplot很適合散佈圖，可以用第三個變數決定顏色，例如我們想呈現SLID這筆資料裡面，200個觀察值的教育與薪資的關係，同時又想觀察跟性別有沒有相關，例如圖4.6：

library(carData)
set.seed(02138)
SLID.sub <- SLID[sample(nrow(SLID), 200),]
qplot(education, wages, data=SLID.sub, size=I(2),shape=I(17),color=sex)

Figure 4.6: qplot散佈圖2

可以看出，性別與教育水準跟薪資之間的關聯並不明顯。
我們也可以加上非線性的迴歸線，圖4.7顯示年齡與薪水呈現倒U的型態，也就是年齡輕與年齡高的薪水相對來得低：

qplot(age, wages, data=SLID.sub, size=I(2),shape=I(4),
      geom=c('point','smooth'))

Figure 4.7: qplot散佈圖3

4.4 qplot多重圖形

圖4.8 顯示四次從SLID的薪資資料的抽樣結果的圖形，疊在同一張圖，方便閱讀。可以看到抽樣的人數越多，分佈都呈現多數人在左邊，越往右邊人數越少：

set.seed(11601)
x1<-c()
    for (i in 1:100){x1[i] <- sample(SLID$wages, 1) 
}
x2<-c()
    for (i in 1:500){x2[i] <- sample(SLID$wages, 1) 
}
x3<-c()
    for (i in 1:1000){x3[i] <- sample(SLID$wages, 1) 
}
x4<-c()
    for (i in 1:2000){x4[i] <- sample(SLID$wages, 1) 
}

mydata <- data.frame(X=c(x1,x2,x3,x4), sample.n=c(rep('100',100),
                rep('500',500),rep('1000',1000),rep('2000',2000)))
mydata$sample.n <- factor(mydata$sample.n, levels=c('100','500',
                                               '1000','2000'))
p<-qplot(X, geom = 'histogram', bins=30, facets = sample.n ~., data=mydata, fill=sample.n)

p + theme_bw()

Figure 4.8: qplot多重圖形

4.4.1 課間練習

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

5 lattice 繪圖

lattice套件用格子(trellis)圖形來產生高品質的資料視覺化，尤其是用在多變量的資料。 Deepayan Sakar在2008年出版「Lattice: Multivariate Data Visualisation in R」，詳細地介紹lattice套件的功能。有興趣的同學可以閱讀。

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

5.1 直方圖 (Bar chart)

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

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

Figure 5.1: lattice直方圖一

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

5.1.1 課間練習

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

5.2 長條圖 (histogram)

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

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

Figure 5.2: lattice長條圖一

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

5.3 雙變數直方圖 (barchart)

lattice並沒有單變數的直方圖，以長條圖代替直方圖。但是比較類別變數與連續變數之間的關係時，可以用直方圖，並且可以控制另一個類別變數。例如圖 5.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 5.3: lattice直方圖二

5.4 表格資料點狀圖 (Dotplot)

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

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

Figure 5.4: lattice點狀圖一

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

如果要畫在同一個座標軸上，需要轉換資料為資料框，如圖 5.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 5.5: lattice點狀圖二

在上面的例子中，

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

5.5 統計資料點狀圖

我們想知道每一類型的樹的平均樹圍，並且畫圖。首先用tapply計算平均值。然後用dotplot畫圓點表示平均樹圍。先用tapply()函數產生每一群樹木的平均值，然後用這個平均值的表格的名稱當作Y軸。請見圖 5.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 5.6: lattice點狀圖一

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

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

Figure 5.7: lattice點狀圖二

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

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

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

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

Figure 5.8: lattice點狀圖三

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

5.6 散佈圖 (xyplot)

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

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

Figure 5.9: lattice散佈圖一

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

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

Figure 5.10: lattice散佈圖二

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

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

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

Figure 5.11: lattice長條圖

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

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

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

Figure 5.12: lattice箱型圖

5.9 折線圖

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

lattice::xyplot(MASS::accdeaths)

Figure 5.13: lattice折線圖一

圖 5.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$。中介變數必須要是類別。

圖 5.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 5.14: 每年的意外死亡折線圖

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

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

6 以ggplot 為基礎的套件

6.1 ggpubr

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

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

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

Figure 6.1: ggpubr長條圖

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

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

Figure 6.2: ggpubr的箱型圖

6.2 ggstatsplot

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

utils::install.packages(pkgs = "ggstatsplot")

我們以ISLR::Default為例，進行交叉列表並卡方檢定，呈現圓餅圖 6.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 6.3: 卡方檢定圓餅圖

兩個連續變數之間可以用相關係數來描述兩者之間的關聯，而ggstatsplot可以顯示相關係數。在執行下面指令前先下載ggside，方法如下：

# downloading the package from GitHub
devtools::install_github(
  repo = "jtlandis/ggside", # package path on GitHub
  dependencies = FALSE,                # assumes you have already installed needed packages
  quick = TRUE                         # skips docs, demos, and vignettes
)

library(ggside)
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 6.4: ggstatsplot散佈圖與迴歸係數

7 函數與ggplot2

還記得一開始的課程目標是畫出類別變數與連續變數之間的關係，並且標示各類別的平均值？請見圖 7.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 7.1: ggplot2散佈圖加上平均值

我們運用stat_summary()計算各類別對應的變數的平均數，還可以加上標準差，並且運用pointrange的指令，顯示平均值與標準差連起來的直線，可以用來比較各類別之間的差異程度。圖 7.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 7.2: ggplot2散佈圖加上平均值及標準差

8 作業

8.1 基礎繪圖指令作業

1. 城市或者州的犯罪率與貧窮程度有關嗎？請針對carData套件裡面的Ericksen這筆資料，製作poverty以及crime的散佈圖，並且加上城市的名稱。

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)。請記得加上圖例。

8.2 ggplot2以及lattice作業

1. 種族與社會經濟地位有關嗎？請用ggplot2繪製「灰階」堆疊圖顯示hsb2的race與ses的關係。

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

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

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是一筆表格資料，有六個系錄取與拒絕男性與女性申請者的人數(Freq)。請轉換這筆資料，然後用畫圖表示六個系男性與女性的錄取率，並回答哪一個系的男性與女性錄取率差異最大。（提示：使用dplyr指令整理資料。圖上面加上資料數值。）

7. 請用lattice畫圖顯示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軸的刻度以及子圖的順序。

9 更新講義時間

最後更新時間: 2024-04-01 19:53:32

社會科學統計方法

資料視覺化

蔡佳泓

2024