本週上課將介紹R
的ggplot2
的套件中的指令,把資料視覺化,雖然之前我們已經看過很多例子,但是仍然有必要好好介紹,包括直方圖、折線圖、散佈圖等等。例如:
library(ggplot2)
students<-read.table('studentsfull.txt', sep='', header=TRUE)
p <- ggplot(students, aes(x=Department, y=Score)) +
geom_point() +
theme_classic() +
stat_summary(fun.y='mean', color='red')
p
有人說,「一張照片勝過千言萬語」。以圖形表現變數之間的關係,在社會科學界經常被學者採用,現在新聞傳播也大量用統計圖表說明,提高讀者的興趣。
R
內建了許多基本的繪圖指令,讓我們可以很輕鬆地使用,但是ggplot2
提升了繪圖的層次,配合之前介紹的dplyr
套件,可以畫出高品質的圖形,所以我們必須學會ggplot2
。另外,lattice
套件提供許多繪圖的功能,也值得介紹。
R
內建了許多基本的繪圖指令,可以很直覺地畫圖描述變數或者是兩個變數之間的關係,使用者同時可以微調一些繪圖的元素。
首先我們先認識一些繪圖的元素。首先是「點」、「線」、「大小」、「粗細」以及「顏色」。
下圖表示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"
## [8] "gray20" "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)
如果要檢視目前繪圖的參數,可以這麼做:
par()
例如我們先設定無邊框、前景為藍色、X, Y軸的座標標記為1.5倍大。在設定之前,先把目前的參數存在一個變數中。
oldpar <- par()
par(bty='n', fg='blue', cex.axis=1.5, col.main='red')
然後我們畫一個直方圖加上常態曲線:
oldpar <- par()
par(bty='n', fg='blue', cex.axis=1.5, col.main='red')
library(foreign)
teds2016<-read.spss('TEDS2016_indQE(CSES).sav', to.data.frame=TRUE)
## re-encoding from CP950
barplot(table(teds2016$Q2)/length(teds2016$Q2), col='brown',
main='Frequencies with Normal Distribution')
curve(dnorm(x, 2, 1), 0, 4, lwd=2, col=3, add=T)
par(oldpar)
在這個例子中,
或者加上一個par(new=TRUE)的指令,但是要注意去掉曲線圖產生的座標以及標題:
barplot(table(teds2016$Q2)/length(teds2016$Q2), col='navy',
main='Frequencies with Normal Distribution')
par(new=T)
curve(dnorm(x, 2, 1), 0, 8, lwd=2, col='orange4',
xaxt='n', yaxt='n', ylim=c(0, 1), ylab='', xlab='')
我們再畫別的資料的變數:
hist(nycflights13::weather$humid)
除了以上常用的參數外,可以參考World is Interesting,裡面介紹更多參數。
我們可以指定R
圖形的顏色,例如red, blue, brown等等。
colors()
另外,有一套16進位制的顏色號碼,可以找到更多顏色(http://www.sthda.com/english/wiki/colors-in-r),以上面的直方圖為例:
par(bg="#CCCCCC", col.lab="#00000F", col.main='#00000F')
barplot(table(teds2016$Q2)/length(teds2016$Q2), col='#FFFF33', main='Distribution of Q2', col.main='#00000F')
curve(dnorm(x, 2, 1), 0, 5, lwd=2, col='#FF3300', add=T)
此外,我們可以使用調色盤的套件以及函數,增添有趣的顏色。例如先下載RcolorBrewer套件:
install.packages("RColorBrewer")
library(RColorBrewer)
barplot(c(2,5,7, 10), col=brewer.pal(n = 4, name = "RdBu"))
實際應用在散佈圖:
with(faithful, plot(waiting, eruptions, col=brewer.pal(n=6, name='RdYlGn')))
那麼我們如何應用到資料分析?例如我們有一筆資料anscombe,其中有兩個連續變數,我們可以用散佈圖來檢視:
with(anscombe, plot(x1, y1, pch=16, cex=2, col='blue',
xlab='X', ylab='Y', main='x1-y1'))
用points有同樣的效果,但是要先幫觀察值畫一個框:
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.260 6.315 7.580 7.501 8.570 10.840
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'))
用points繪圖標示兩個變數對應的數值,可以控制散佈圖顯示的圖案:
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=16, cex=2, col='darkblue'))
abline(v=8, lwd=2, lty=3)
上圖顯示x1超過8是深藍色的圓點,小於8是紅色的點。
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))
用text可以加上文字,讓讀者更容易瞭解觀察值的屬性。pos 控制文字的位置。也可以貼上觀察值本身的數字:
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))
繪圖時要記得兩個變數的長度需要一致。請大家幫忙想想如何幫上面這個圖裡面的觀察值數字前後加上括號,中間要有逗號?
直方圖可以表現單一類別變數的分佈,例如:
students<-read.table('studentsfull.txt', header=T, sep='')
stu <- table(students$Department)
barplot(stu, main="Departments",
xlab="", ylab="frequency")
直方圖的對象應該是具有名稱的摘要,所以要先把數字、字串或者是類別的向量改為表格的形式。直方圖也可以改直方本身的顏色以及外框的顏色。
barplot(100*stu/nrow(students), main="Departments",
xlab="", ylab="Percent", border='red', col='cyan', ylim=c(0, 20))
如果想要表現兩個變數的交集,例如經濟系裡面有多少男生或女生的比例,可以先產生一個交叉列表,然後用prop.table()的函數產生條件機率,就可以繪出堆疊的直方圖:
student.table <- table(students$Gender, students$Department)
barplot(100*prop.table(student.table, margin=2),
col=c('brown', 'white'),
legend = levels(unique(students$Gender)))
如果想要更改類別的名稱,可以直接在barplot()之中設定參數:
barplot(100*prop.table(student.table, margin=2),
names.arg=c("Aer","Che.", "Eco.", "Eng.", "Jou.", "Mec.","Phy."))
箱形圖於1977年由普林斯頓大學統計系教授約翰·圖基(John Tukey)發明。可以表現單一連續變數的中位數、25分位數、75分位數、極端值等等,也可以比較不同類別之下,連續變數的分佈情形。
例如我們觀察mpg這個變數的盒型圖:
boxplot(mtcars$mpg, ylim=c(0,40), yaxt='n')
中間的粗線代表中位數。盒型圖上方的線稱為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.275 19.200 22.800
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.425 19.200 22.800
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.8625
## / lower 4.3625
# upper inner fence
min(max(mtcars$mpg), upper)
## [1] 33.8625
# lower inner fence
max(min(mtcars$mpg), lower)
## [1] 10.4
畫箱形圖對照上述的計算結果:
boxplot(mtcars$mpg, ylim=c(0,40), yaxt='n')
axis(2, at=c(1:40, by=5), labels=c(1:40,by=5))
上述的資料並沒有極端值,我們可以觀察美國各州的面積,其中有幾個州特別的大:
boxplot(state.area, ylab="Area of State")
然後請輸入下列指令:
identify(rep(1, length(y)), y, labels=seq_along(y)))
用滑鼠對著圖形點擊,如果出現已經找到最近的點,可以按旁邊的finish,然後會呈現哪幾個觀察值被標示在圖形上面。例如第二個州是阿拉斯加,面積最大。
我們可以對類別變數繪製另一個連續變數的箱型圖,例如:
Orange$tree <-ordered(Orange$Tree, levels=c(1,2,3,4,5))
with(Orange, boxplot(circumference ~ tree))
由上圖可以看出,第一類型的四分位距比較小,其次是第三類型。而第四類型的中位數最大。
因為連續變數有一定數量的值,而非像類別變數的值是離散的,所以應該用長條圖來呈現。例如:
hist(USArrests$Assault, density=6)
上圖的指令中,density參數可以產生斜線陰影。而連續變數的每一個值的相對比例可以幫助我們了解分佈型態,因此,freq=F參數強制長條圖呈現相對的比例,而非次數。
hist(USArrests$Assault, col="red", freq = F, xlab="Assault", main="Assault in 50 States")
我們用模擬的資料來表現長條圖的參數break的用途,左邊的圖最多有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))
在長條圖上面可以加上特定分布的曲線,例如常態分佈曲線或者是一致分佈曲線:
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)
從上面的圖可以看出模擬的資料近似常態分佈,而常態分佈曲線比一致分佈曲線更接近資料分佈。
之前在介紹散佈圖時,曾經設定參數type=‘n’,使觀察值不顯示出來,而除了type=‘n’,還有b, c, h, o, p等5個選項。而且,可以用單一變數繪圖表示分佈型態。例如:
with(airquality, plot(Wind, type='b'))
或者是:
plot(state.area, type='o', pch=16, cex=1.2, lty=2, lwd=2, col='red')
可以設定線條的形式、粗細以及連接點的顏色、大小與形狀。
資料有時間序列的型態,適用以折線圖顯示,例如LakeHuron這筆資料是時間序列資料,裡面的參數設定為type=‘o’。
plot(LakeHuron, type = "o", pch=16, cex=1.2, lty=2) ## Index plot
那麼我們要如何建立時間序列資料?ts()函數可以轉換一個資料框為時間序列資料,要注意設定起始的時間點,可以以月份來切分,如果頻率為1,就是以1年為單位,如果頻率為2,就是半年,頻率為4,就是3個月:
trend<-read.csv("tondutrend.csv",
header=T, sep=",",fileEncoding="BIG5")
tonduts<-ts(trend, start=1992.6,frequency=2)
相對於其他軟體,R
的優點在於允許時間點有缺漏,畫出來的圖形就有缺少某個或是某些時間點,不會強制連接前後兩個時間點。接下來是畫折線圖:
# Plot
par(xpd=NA, mar=par()$mar+c(2.5, 0, 0, 0))
plot(tonduts, plot.type=c("single"), lty=c(1,2,2,3),
ylab="%",xlab=NULL,pch='1', lwd=3,frame.plot=F,
col=c("gray20","gray60", "black", "gray80"),
xaxt="n")
axis(1, at=seq(1992,2014,by=2))
axis(2, at=seq(10,60,by=10))
legend("bottomright", c("unification","status quo"),
inset=c(0.35, -0.4), col=c("gray20","gray60"),
lty=c(1,2), bty='n', lwd=3)
legend("bottomright", c("independence","don't know"),
inset=c(0, -0.4), col=c("black", "gray80"),
lty=c(2,3), bty='n', lwd=3)
text(2000, 30, paste("First Party Turnover"))
散佈圖可表示單變數的分佈,也可以呈現兩個變數之間的關係。例如rivers這筆資料是美國河流的長度,我們由大排到小排序變數之後,畫成散佈圖:
plot(rev(sort(rivers)), type='l', col='darkblue', lwd=1.5)
LakeHuron這筆資料是這個湖的水位變化的時間序列資料,也可用散佈圖表示:
plot(LakeHuron, col='darkgreen', lwd=1.5)
兩個變數的散佈圖可以這樣畫:
plot(Murder ~ Assault, data=USArrests, cex=1.5)
散佈圖中的圓點可變成特定的字串:
plot(Murder ~ Assault, data=USArrests, pch='x', cex=1.5)
散佈圖中的圓點可變成向量的資料,例如:
plot(Murder ~ Assault, data=USArrests, pch=row.names(USArrests))
圖上的點變成一個字母,代表各州州名的第一個字母。 如果想加上完整的字串,可用text()函數:
plot(Murder ~ Assault, data=USArrests, type='n')
with(USArrests, text(Assault, Murder, labels=row.names(USArrests), col='#CC3366'))
散佈圖可加上迴歸線:
model <- lm(Murder ~ Assault, data=USArrests)
plot(Murder ~ Assault, data=USArrests, cex=1)
abline(model, col='red', lty=2, lwd=1.5)
例如我們先排序某一個變數,然後畫圓點圖:
dt <- UsingR::carsafety[order(UsingR::carsafety$Driver.deaths),] # sort
par(mfrow=c(1,2))
dotchart(dt$Driver.deaths, main='Dotchart')
plot(dt$Driver.deaths, main='Plot')
接下來我們嘗試在圖形中加入分組排序:
dt <- UsingR::carsafety[order(UsingR::carsafety$Driver.deaths),] # sort
dt$type <- factor(dt$type) # it must be a factor
dt$color[dt$type=='large'] <- "red"
dt$color[dt$type=='minivan'] <- "blue"
dt$color[dt$type=='midsize'] <- "darkgreen"
dt$color[dt$type=='subcompact'] <- "black"
dt$color[dt$type=='compact'] <- "orange"
dt$color[dt$type=='SUV'] <- "brown"
dotchart(dt$Driver.deaths,labels=row.names(dt),cex=.7,groups= dt$type,
main="Number of driver death\ngrouped by Type",
xlab="Number", gcolor="black", color=dt$color)
symbols()可以產生指定形狀、大小的散佈圖,而且可以進一步根據另一個變數的觀察值,調整散佈點的大小,例如:
with(anscombe, symbols(x2, y1, circles=anscombe$x1,
inches=0.2, fg='blue'))
可以看出有些圈圈比較大,有些比較小。再舉一個例子:
with(USArrests, symbols(Murder, Assault, circles=UrbanPop,
inches=0.12, bg="red"))
上圖顯示,雖然搶案與攻擊成正比,但是都市化人口越多的地方,不見得有更多的類似案件。
接下來我們介紹用ggplot2
繪圖。基本的指令為:
ggplot(data, aes(x, y, group,…)) + geom_object () + theme()
ggplot2
的例子:http://r-statistics.co/Top50-Ggplot2-Visualizations-MasterList-R-Code.html#Ordered%20Bar%20Chart
變數屬於數字、類別或是字串都可以以直方圖表示。例如有一筆資料是雞的成長數字,先看一下資料:
head(ChickWeight)
## weight Time Chick Diet
## 1 42 0 1 1
## 2 51 2 1 1
## 3 59 4 1 1
## 4 64 6 1 1
## 5 76 8 1 1
## 6 93 10 1 1
如果我們想畫雞的飼料類型的分佈,可以這樣畫:
ggplot(ChickWeight, aes(x=Diet), stat='bin') +
geom_bar(fill='red', col='blue') +
ggtitle('Diet')
注意上面的設定為stat=“bin”,代表計算每一個類別發生的次數。如果資料是統計過的,也就是表格,要設定stat=“identity”。首先我們轉換資料為表格,再轉換為資料框:
CW <-table(ChickWeight$Diet)
CW1 <- as.data.frame.table(CW)
CW1
## Var1 Freq
## 1 1 220
## 2 2 120
## 3 3 120
## 4 4 118
然後畫直方圖:
ggplot(CW1, aes(x=Var1, y=Freq)) +
geom_bar(fill='navy', stat='identity') +
ggtitle('Diet')
又例如統計學生資料中各系所的人數:
students<-read.table('studentsfull.txt', 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='darkgreen') +
ylab('Count')
g1
Table and Barplot
加上數據
用已經統計好的數字當作標記(label),而該數字加上或是減去某個單位當作y,用geom_text()標上數據:
g1 + geom_text(data=stu.t, aes(x=Dep, label=Freq, y=Freq+0.2), size=8)
之前已經畫過有兩個變數的直方圖,我們用mutate產生一個新的變數Pass,兩者的交叉列表可以顯示如下:
library(dplyr)
students<-read.table('studentsfull.txt', header=T, sep='')
stu <- dplyr::mutate(students, Pass=(Score>70))
head(stu)
## ID Name Department Score Gender Pass
## 1 10322011 Ariel Aerospace 78 F TRUE
## 2 10325023 Becky Physics 86 F TRUE
## 3 10430101 Carl Journalism 69 M FALSE
## 4 10401032 Dimitri English 83 M TRUE
## 5 10307120 Enrique Chemistry 80 M TRUE
## 6 10207005 Fernando Chemistry 66 M FALSE
然後畫成並列的直方圖:
ggplot(stu, aes(x=Department, fill=Pass)) +
geom_bar(position='dodge')
Barplot Side by Side
之前的直方圖的色塊表示條件機率,我們可以用summarize產生Count這個變數,代表每個系通過或是沒通過的個案數。然後用mutate產生以系所為邊際機率的條件機率。
stu <- dplyr::mutate(students, Pass=(Score>70))
stu.ag <- summarize(group_by(stu, Department, Pass), Count=n())
stu.ag <- mutate(stu.ag, Pct=Count/sum(Count))
head(stu.ag)
## # A tibble: 6 x 4
## # Groups: Department [4]
## Department Pass Count Pct
## <fct> <lgl> <int> <dbl>
## 1 Aerospace TRUE 4 1
## 2 Chemistry FALSE 1 0.333
## 3 Chemistry TRUE 2 0.667
## 4 Economics FALSE 2 0.5
## 5 Economics TRUE 2 0.5
## 6 English TRUE 3 1
畫成圖表示上列數據:
ggplot(stu.ag, aes(x=Department, y=Pct, fill=Pass)) +
geom_bar(stat='identity') +
scale_y_continuous(label = scales::percent)
Barplot with Percentage
ggplot2
的長條圖跟直方圖類似,但是適用在連續變數。參數stat=‘bin’指定每一長條的寬度,寬度越寬,越容易黏合在一起,但是也可能越不容易分辨分佈的型態:
ggplot(mtcars, aes(x=mpg)) + geom_bar(stat='bin', binwidth=1)
Histogram
長條圖可以轉換成密度圖,而且可以按照對應的類別變數並列在同一張圖:
library(ggplot2)
ggplot(mtcars, aes(x=wt, fill=as.factor(am))) +
geom_density(position="identity", alpha=.4)
Density Histogram
如果需要表現兩個變數之間的關係,散佈圖可以顯示兩個變數對應的位置:
ggplot(mtcars, aes(x=wt, y=mpg)) +
geom_point()
Scatter Plot
可以加上迴歸線表現兩個變數的相關:
ggplot(mtcars, aes(x=wt, y=mpg)) +
geom_point(col='red') +
geom_smooth(method="lm", se=F, col='blue')
Scatter Plot with Regression Line
也可以用非線性迴歸線表示兩者的相關:
ggplot(Orange, aes(x=age, y=circumference)) +
geom_point(aes(col=Tree), size=3) +
geom_smooth(method="loess", se=F) +
theme_bw()
Scatter Plot with LOESS Line
兩個變數的散佈圖可表示兩個變數之間的相關:
ggplot(anscombe, aes(x=x1, y=y1)) +
geom_point(size=3)
Scatter plot
如果加上另一個變數,可以表示當另一個變數等於某一個類別時,兩個變數是否仍然維持一樣的關係?我們用shape設定形狀隨著類別變數而改變,然後隨著類別變數改變顏色:
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)
Scatter plot with a control variable
也可以隨著連續變數而改變顏色深淺:
sc2<-ggplot(anscombe, aes(x=x1, y=y1))
sc2 + geom_point(aes(color=x2) , size=3)
Scatter plot with a contineous variable
盒型圖可以表現一個變數或者是一個類別變數對應另一個連續變數的分佈,不過ggplot2
無法直接畫出一個連續變數的盒型圖,所以我們先創造一個只有一個值的變數,然後對應我們要顯示的連續變數,例如:
mtcars <- mutate(mtcars, X=1)
ggplot(mtcars, aes(x=X, y=mpg)) +
geom_boxplot() +
labs(x="",y='mpg') +
stat_summary(fun.y=mean, geom="point", shape=16, size=1.5, col='red2') +
theme_bw()
在上面這個例子,fun.y=mean可以產生平均值,與中間的中位數不同。
盒型圖可以用來比較一個類別變數對應的連續變數的分佈,例如我們想要顯示不同的氣缸數(cyl)的馬力(hp),可以把y設定為(hp),然後把x, group設定為(cyl):
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
Boxplot of Two Variables
由於(cyl)是數字變數,我們發現有一些我們不需要的類別,沒有資料與其對應,因此我們用由於(scale_x_discrete)加以調整,也就是只留下需要的類別數字,然後加上所需要的標籤:
g1 + scale_x_discrete(limit = c(4, 6, 8),
labels = c("4","6","8")) +
theme_bw()
類似的技巧可以適用在Y軸,例如我們只想顯示50到200的馬力,超過200以上的資料就會被去掉:
g1 + scale_y_continuous(limits = c(50,200)) +
scale_x_discrete(limit = c(4, 6, 8),
labels = c("4","6","8")) +
theme_bw()
過去有所有的圓餅圖,甜甜圈圖可以顯示資料的次數分佈,例如我們想知道ggplot2
裡面的的diamonds
資料中的cut的次數分配,先針對這個變數計算個數:
#
dt <-group_by(diamonds, cut)
dat = summarize(dt, count=n())
dat
## # A tibble: 5 x 2
## cut count
## <ord> <int>
## 1 Fair 1610
## 2 Good 4906
## 3 Very Good 12082
## 4 Premium 13791
## 5 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))
dat
## # A tibble: 5 x 5
## cut count fraction ymax ymin
## <ord> <int> <dbl> <dbl> <dbl>
## 1 Fair 1610 0.0298 0.0298 0
## 2 Good 4906 0.0910 0.121 0.0298
## 3 Very Good 12082 0.224 0.345 0.121
## 4 Premium 13791 0.256 0.600 0.345
## 5 Ideal 21551 0.400 1 0.600
然後畫圖:
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
Ring Plot
之前我們用R
的基本指令對airquality
這筆時間序列資料畫折線圖,
ggplot(airquality, aes(x=Day, y=Wind)) +
geom_line()
Line Plot with Wrong Index
這個圖看起有很多重疊的地方,原因是日期是從1日到30日或是31日循環出現,所以我們可以創造一個新的變數,可以顯示時間序列:
N <- summarize(airquality, n())
airquality <- mutate(airquality, index=seq_along(1:N[,1]))
ggplot(airquality, aes(x=index, y=Wind)) +
geom_line()
Line Plot with an Index
我們示範如何把原來的月份與日變數轉換成日期:
airquality$Date <- paste0(airquality$Month, "/", airquality$Day)
airquality$d<-as.character(airquality$Date)
airquality$d<-as.Date(airquality$d, format='%m/%d')
然後用日期畫折線圖,X軸標示的5, 6,…,10為月份:
ggplot(airquality, aes(x=d, y=Wind)) +
geom_line() +
labs(x='Month')
Line Plot with Date
進一步可以顯示不只一項的折線圖,但是我們需要轉置資料為長表,在轉置之前先選擇兩個變數加上日期:
airquality.n <- select(airquality, Date, Wind, Temp)
library(reshape2)
airquality.l<-melt(airquality.n, id.vars=c('Date'))
由於這兩個變數的大小範圍接近,可以畫在一張圖上面:
ggplot(airquality.l, aes(x=Date)) +
geom_line(aes(y=value, col=variable))
## geom_path: Each group consists of only one observation. Do you need to
## adjust the group aesthetic?
Multiple Timeseries in a Line Plot
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.4221 -1.7924 -0.3788 1.2249 5.5317
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.002875 2.642659 12.867 2.82e-13 ***
## wt -2.878575 0.904971 -3.181 0.003574 **
## hp -0.037479 0.009605 -3.902 0.000546 ***
## am 2.083710 1.376420 1.514 0.141268
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.538 on 28 degrees of freedom
## Multiple R-squared: 0.8399, Adjusted R-squared: 0.8227
## F-statistic: 48.96 on 3 and 28 DF, p-value: 2.908e-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分別產生線段:
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
然後加上代表係數的圓點,並加上係數的名稱,以及把X軸與Y軸對調,就大功告成。
s1 + geom_point(data=dt, aes(x=row.names(dt), y=coef), size=4, shape=21, fill="white") +
scale_x_discrete(breaks=c("1","2","3"),
labels=c("Weight", "Horse Power", "Auto")) +
ylab("") +
xlab("Estimates") +
coord_flip() +
theme(axis.text.y = element_text(face="bold", color="#993333",
size=14) )
Coefficient Plot
從這個係數圖可以看出是否自動排檔對於依變數的影響相對來得大,但是有可能等於0,因為標準誤的下限低於0。
雖然上圖已經相當美觀,但是仔細觀察,會發現係數的順序反過來了。要改變這個情形,可以用rev()這個函數改變向量的元素順序,例如:
V<-c(50, 100, 30, 200)
rev(V)
## [1] 200 30 100 50
除了上述的圖形,ggplot2
還有更多有趣的圖形,也允許使用者設定圖形風格、資料點顏色、副標題、圖例的標題等等,請同學上網多多參考他人的例子。例如以下的圖顯示經濟學人雜誌風格的直方圖:
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
Deepayan Sakar出版「Lattice: Multivariate Data Visualisation in R」,詳細地介紹lattice套件的功能。有興趣的同學可以自行上網下載。
我們以Figure 1顯示星際大戰中的演員膚色的直方圖:
library(lattice); library(dplyr)
barchart(starwars$skin_color,col="gray60",
scale=list(cex=0.9))
長條圖可用來描述中央趨勢:
histogram(faithful$eruptions, col='gray60')
histogram跟R
基礎指令hist函數得到的圖形相比,多了一個外框。
lattice並沒有單變數的直方圖,以長條圖代替直方圖。但是比較類別變數與連續變數之間的關係時,可以用直方圖,並且可以控制另一個類別變數。例如:
library(MASS)
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
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=brewer.pal(n=4, 'RdGy'),
main=paste0("Crab's Body Depth(mm)"))
我們可以直接對一個表格資料畫圖,表示兩個變數之間的關係。例如使用VADeaths這筆表格資料,畫圖如下:
lattice::dotplot(VADeaths, group=F, type="o")
以上的Figure 2顯示,四群觀察值的年齡對應的死亡率不太一樣。
如果要畫在同一個座標軸上,需要轉換資料為資料框:
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))
在上面的例子中,
兩個連續變數可以用散佈圖表示相對位置:
with(lattice::ethanol, xyplot(NOx ~ E))
用途與plot()一樣,只是內建的顏色不同。但是xyplot()可控制另一個變數:
with(lattice::ethanol, xyplot(NOx ~ E|C, col='#FF3300'))
在lattice,類別變數對應連續變數的密度長條圖,可以一次展示在同一個圖上面。用以下指令畫圖:
histogram(~temp|origin, data=nycflights13::weather, type="density", layout=c(1,3))
箱形圖也可以表示連續變數在每一個類別上面的分佈:
bwplot(as.factor(cyl) ~ mpg, data=mtcars)
以上介紹的lattice套件繪圖功能,雖然品質不像ggplot2,但是提供R
以外的選擇,如果經常運用,會有不錯的效果。
mtcars
這筆資料的wt以及mpg繪製散佈圖,並且加上車輛的名稱:
ChickWeight
的,Diet變數畫出直方圖,並且標示次數(提示:text(x, y, labels),x, y 分別是標示次數的位置,而標示次數要轉為字串)。
mtcars
這筆資料中的變速箱(am)是屬於自動或是手動排檔之中,gear的比例。
mtcars
這筆資料中的gear的各個類別的mpg。
ggplot2
繪製直方圖顯示hsb2的race與ses的關係。
flights
這筆資料之中的arr_delay,畫圖顯示每一家航空公司(carrier)的分布情況:
airquality
這筆資料的Ozone
變數折線圖:
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))
DT
## # A tibble: 5 x 4
## Country Populist Mixed Nonpopulist
## <chr> <dbl> <dbl> <dbl>
## 1 Spain 52 70 85
## 2 France 58 72 84
## 3 Germany 65 70 88
## 4 Italy 50 65 68
## 5 Sweden 73 81 90
ggthemes
套件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)