1 課程目標

本週上課將介紹Rggplot2的套件中的指令,把資料視覺化,雖然之前我們已經看過很多例子,但是仍然有必要好好介紹,包括直方圖、折線圖、散佈圖等等。例如:

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


2 基本指令

有人說,「一張照片勝過千言萬語」。以圖形表現變數之間的關係,在社會科學界經常被學者採用,現在新聞傳播也大量用統計圖表說明,提高讀者的興趣。
R內建了許多基本的繪圖指令,讓我們可以很輕鬆地使用,但是ggplot2提升了繪圖的層次,配合之前介紹的dplyr套件,可以畫出高品質的圖形,所以我們必須學會ggplot2。另外,lattice套件提供許多繪圖的功能,也值得介紹。

2.1 基本概念

R內建了許多基本的繪圖指令,可以很直覺地畫圖描述變數或者是兩個變數之間的關係,使用者同時可以微調一些繪圖的元素。

首先我們先認識一些繪圖的元素。首先是「點」、「線」、「大小」、「粗細」以及「顏色」。

  • cex 控制點的大小
  • col 控制顏色,也可以是灰階的濃度
  • pch 控制點的形狀
  • lwd 控制線的粗細
  • lty 控制線的形狀
  • 下圖表示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)

    從以上的指令以及圖形不難看出一些常用圖形參數。另外要說明幾個參數:
  • 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軸標題的大小
  • cex.axis控制X, Y軸刻度的大小
  • cex.name控制直方圖的X或 Y軸的類別文字的大小
  • bg 設置繪圖區背景色,例如bg=‘red’或是’blue’, ’gray20’等等。
  • fg 設置繪圖區前景色,包括點、線以及刻度顏色。
  • bty 設置繪圖邊框的形式。= “o”(內建):四周邊框;= “l”:左下邊框;= “7”:右上邊框;= “c”:上左 下邊框;= “u”:左下右邊框;= “]”:上右下邊框;= “n”:無邊框。
  • 以上這些參數可以在每一次繪圖時設定,也可以一開始就用par()函數設定,讓每一次繪圖都會得到相同的背景顏色、文字顏色、大小等等。

    如果要檢視目前繪圖的參數,可以這麼做:

    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)

    在這個例子中,

    • 先計算Q2這個變數每一類別的百分比,然後畫成直方圖
    • 加上常態曲線,其中平均值為2,標準差為1

    或者加上一個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,裡面介紹更多參數。

    2.2 顏色

    我們可以指定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))

    繪圖時要記得兩個變數的長度需要一致。請大家幫忙想想如何幫上面這個圖裡面的觀察值數字前後加上括號,中間要有逗號?


    2.3 直方圖 (Bar Plot)

    直方圖可以表現單一類別變數的分佈,例如:

    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.")) 

    2.4 箱形圖 (Box Plot)

    箱形圖於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")

    然後請輸入下列指令:

      y<-state.abb

    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))

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


    2.5 長條圖 (Histogram)

    因為連續變數有一定數量的值,而非像類別變數的值是離散的,所以應該用長條圖來呈現。例如:

    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)

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


    2.6 折線圖 (Line Chart)

    之前在介紹散佈圖時,曾經設定參數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')

    可以設定線條的形式、粗細以及連接點的顏色、大小與形狀。

    2.7 時間序列資料

    資料有時間序列的型態,適用以折線圖顯示,例如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"))

    在這個圖形中,我們設定xpd=NA,允許圖形超過界線,而且設定圖形的邊線可以到原來區域的底部2.5個文字。這是因為折線有四條,需要一定的空間容納圖例,而為了避免與圖形重疊起見,我們把圖例放在X軸右下方,所以除了要設定參數“bottomright”,還要設定inset(0.35, -0.6)以及inset(0, -0.6)。除了設定圖例出現的位置,還可以設定bty=“n”,強制圖例沒有外框。
    最後,我們加上一點註解文字到圖形上面,當然,我們可以加上線條或是箭頭,請大家研究一下segment這個指令。
    相信這個折線圖會讓你的報告比其他軟體畫的圖看起來更專業。

    2.8 散佈圖

    散佈圖可表示單變數的分佈,也可以呈現兩個變數之間的關係。例如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)

    在這個例子中:

    • 先建立一個迴歸模型,然後用abline()把迴歸線加入散佈圖。
    • 注意迴歸模型的Y在左邊,X在右邊,跟plot()裡面的順序一樣。

      2.8 圓點圖(dotchart)

      圓點圖類似散佈圖,以圓點表示資料的分佈。但是圓點圖會加上虛線,引導視覺到座標軸。而且圓點圖把變數的值對應在X軸,散佈圖則在Y軸。

    例如我們先排序某一個變數,然後畫圓點圖:

    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)

    在這個例子中:

    • 針對駕駛死亡人數進行排序
    • 建立六種車款的顏色
    • 根據六種車款畫圓點圖表示每一品牌汽車的駕駛死亡率

    2.9 特殊點狀圖

    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"))

    上圖顯示,雖然搶案與攻擊成正比,但是都市化人口越多的地方,不見得有更多的類似案件。


    3 圖形進階(ggplot2)

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

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

  • geom_object可以是geom_bargeom_pointgeom_rect等等。
  • 繪圖的對象應該是資料框或是tibble的變數之一。
  • 繪圖之前應該特別注意變數是哪一種型態。
  • 這個網頁有許多ggplot2的例子:http://r-statistics.co/Top50-Ggplot2-Visualizations-MasterList-R-Code.html#Ordered%20Bar%20Chart

  • 3.1 直方圖 (Bar plot)

    變數屬於數字、類別或是字串都可以以直方圖表示。例如有一筆資料是雞的成長數字,先看一下資料:

    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

    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

    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

    Barplot with Percentage


    3.2 長條圖 (Histogram)

    ggplot2的長條圖跟直方圖類似,但是適用在連續變數。參數stat=‘bin’指定每一長條的寬度,寬度越寬,越容易黏合在一起,但是也可能越不容易分辨分佈的型態:

    ggplot(mtcars, aes(x=mpg)) + geom_bar(stat='bin', binwidth=1)
    Histogram

    Histogram

    長條圖可以轉換成密度圖,而且可以按照對應的類別變數並列在同一張圖:

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

    Density Histogram


    3.3 散佈圖 (scatter plot)

    如果需要表現兩個變數之間的關係,散佈圖可以顯示兩個變數對應的位置:

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

    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

    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

    Scatter Plot with LOESS Line


    3.4 三個變數之散佈圖

    兩個變數的散佈圖可表示兩個變數之間的相關:

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

    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

    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

    Scatter plot with a contineous variable


    3.5 盒型圖 (Box plot)

    盒型圖可以表現一個變數或者是一個類別變數對應另一個連續變數的分佈,不過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

    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()


    3.6 甜甜圈圖

    過去有所有的圓餅圖,甜甜圈圖可以顯示資料的次數分佈,例如我們想知道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

    Ring Plot


    3.7 折線圖

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

    ggplot(airquality, aes(x=Day, y=Wind)) +
      geom_line()
    Line Plot with Wrong Index

    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

    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

    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

    Multiple Timeseries in a Line Plot

    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.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這個參數,產生線段,原理是xxend,以及yyend分別產生線段:

    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

    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


    4 lattice 套件

    Deepayan Sakar出版「Lattice: Multivariate Data Visualisation in R」,詳細地介紹lattice套件的功能。有興趣的同學可以自行上網下載。

    4.1 直方圖 (Bar chart)

    我們以Figure 1顯示星際大戰中的演員膚色的直方圖:

    library(lattice); library(dplyr)
    barchart(starwars$skin_color,col="gray60",
                scale=list(cex=0.9))

    4.2 長條圖 (histogram)

    長條圖可用來描述中央趨勢:

    histogram(faithful$eruptions, col='gray60')

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

    4.3 直方圖 (barchart)

    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)"))

    4.4 點狀圖 (Dotplot)

    我們可以直接對一個表格資料畫圖,表示兩個變數之間的關係。例如使用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))

    在上面的例子中,

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

    4.5 散佈圖 (xyplot)

    兩個連續變數可以用散佈圖表示相對位置:

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

    用途與plot()一樣,只是內建的顏色不同。但是xyplot()可控制另一個變數:

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

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

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

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

    4.7 類別變數與連續變數的箱形圖

    箱形圖也可以表示連續變數在每一個類別上面的分佈:

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

    以上介紹的lattice套件繪圖功能,雖然品質不像ggplot2,但是提供R以外的選擇,如果經常運用,會有不錯的效果。


    5 作業

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

    2. 請針對ChickWeight的,Diet變數畫出直方圖,並且標示次數(提示:text(x, y, labels),x, y 分別是標示次數的位置,而標示次數要轉為字串)。

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

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

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

    6. 請畫圖表示mtcars這筆資料中的gear的各個類別的mpg

    7. 請用USArrests這筆資料,畫圖表示每一州的謀殺與攻擊之間的關係,並且加上州的縮寫。

    8. 種族與社會經濟地位有關嗎?請用ggplot2繪製直方圖顯示hsb2raceses的關係。

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

    10. 請畫出airquality這筆資料的Ozone變數折線圖:

    11. 請畫點狀圖顯示下列的資料(原始出處與圖形請見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))
    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

    提示:

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

    12. 請畫出以下四筆資料的密度長條圖在同一張圖上(提示:運用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)

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

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