關卡 1

這個課程的目是讓同學理解,如何依據數據的型態,在一個維度上適當的呈現數據。

關卡 2

請同學先打開plot的說明。

?plot

關卡 3

plot主要有兩個參數:xy。同學可能會認為,x代表x座標的值,y代表y座標的值。一般來說這是對的,但是也不盡然。plot是個會依據x的型態不同、y的型態不同,而有不同行為的函數。

關卡 4

首先,我們拿R的內建的infert資料集為範例。這是當年研究流產與不孕症之間關係的資料集。請同學輸入:?infert先看看資料對每個欄位的說明。

?infert

關卡 5

infert得第一欄(education)代表該名婦女的教育程度;第二欄(age)代表年齡;第三欄(parity)代表已經生產的子女數(經產狀況);第四欄(induced)代表人工流產的次數;第五欄(case)代表是實驗組或是對照組(也就是有無不孕);第六欄(spontaneous)代表自然流產的次數;第七欄與第八欄則是與實驗設計相關的資料。

關卡 6

接著我們來具體看看資料。請同學輸入:View(infert)來看看實際的資料的值。這裡每筆資料都代表一名婦女。

View(infert)

關卡 7

我們等等先拿infert$spontaneous做說明。請同學先把infert$spontaneous存到spon變數中。

spon <- infert$spontaneous

關卡 8

請問同學,spon是什麼型態?

class(spon)

關卡 9

在R中,plot是最泛用的繪圖指令。我們先請同學輸入:plot(spon)。

plot(spon)

關卡 10

plot收到一個數值向量或整數向量時,會把該向量的值當成y軸座標,而對應的x軸座標則為數據的順序。舉例來說,這張圖上的第一筆數據的座標是:(1,spon[1])、第二筆數據的座標是:(2,spon[2])以此類推。同學可以搭配之前打開的infert資料來看。最左邊的點,x軸座標是1代表是第一筆婦女的資料,y軸座標是2.0,代表她自然流產過兩次。

關卡 11

這張plot(spon)的結果,可以讓我們觀察spon是否和資料的順序有相關性。由於整個圖形並沒有看出特別的趨勢,所以「沒有證據」顯示兩者有特別的關係。請大家記下這張圖的特徵:數據均勻的上下散布,並沒有群聚的現象。這就是「沒有證據」的徵狀。

關卡 12

接著,我們請同學先輸入:spon<-factor(spon),轉換一下數據形態。

spon <- factor(spon)

關卡 13

我們再輸入一次:plot(spon)

plot(spon)

關卡 14

plot收到一個類別向量或字串向量時,會依序畫出不同類別的個數。因為這種狀況下,R不知道合理的方法把數據轉換成y座標(數值)。舉例來說,圖中從左邊數來第一個灰色的bar,x軸的說明文字是:0,y軸座標約140,代表spon中有140個0

關卡 15

請問同學,spon中大約有多少個1?

70

關卡 16

請同學輸入:table(spon)。這個函數會計算一個向量中各種值的個數。

table(spon)

關卡 17

R會告訴各位同學,spon中剛好有71個1

關卡 18

我們也拿spon來複習一下初等統計學的數值系統。根據?infert的說明,spon的值為0時代表該名婦女沒有自然流產過,值為1時代表該名婦女自然流產過一次,值為2時代表該名婦女自然流產過兩次以上。請問同學,spon比較適合名目尺度(選A)、順序尺度(選B)、區間尺度(選C)、還是比例尺度(選D)?

B

關卡 19

以自然流產數據spon為例,我也認為第二張的長條圖(把數據當成順序尺度),是比較容易解讀的。事實上,統計學家在看名目尺度與順序尺度的資料時,也只能看各種類別的個數。而順序尺度應該在成像長條的排列順序時,依照類別的順序做排序。

關卡 20

這張圖直接告訴我們:spon只有3種值,而「0次流產」的婦女數量比「1次流產」和「超過2次流產」的次數還要多。

關卡 21

另外一種常用的圖形是圓餅圖。請同學輸入:pie(table(spon))

pie(table(spon))

關卡 22

圓餅圖利用面積與角度,告訴我們數據分佈的比率。由圖中可以明顯看出,“0”類別的個數佔了整體數據的一半。這點是長條圖比較難直接看出來的。R的pie函數是需要名字與個數比率的,所以我們需要先利用table計算每個類別的個數,再把table的結果傳遞給pie使用。

關卡 23

請同學把infert$age存到變數age

age <- infert$age

關卡 24

請問age的型態是?

class(age)

關卡 25

請同學輸入plot(age)

plot(age)

關卡 26

這張plot(age)的結果,可以讓我們觀察age是否和資料的順序有相關性。由於整個圖形並沒有看出特別的趨勢,所以「沒有證據」顯示兩者有特別的關係。請大家記下這張圖的特徵:數據均勻的上下散布,並沒有群聚的現象。這就是「沒有證據」的徵狀。

關卡 27

我們除了畫散佈圖,也可以畫折線圖。請同學輸入:plot(age,type="l")

plot(age, type = "l")

關卡 28

透過折線圖,我們可以看到整體數據有個週期性。整個數據有一個不甚明顯的三個波峰。如果我們繼續鑽研這個數據,會發現數據中其實有許多組比對組合,我推測這三個不明顯的趨勢,可能和資料的收集有關。

關卡 29

分析時,我們常常用直方圖來觀看區間尺度或比例尺度的數據分佈。請同學輸入:x<-hist(age)來繪製age的直方圖,觀察age的分佈。

x <- hist(age)

關卡 30

根據右邊的直方圖,我們看到有三個年齡層的樣本數特別高。雖然圖中可以告訴我們,最高峰位於25~30歲之間,但是並沒辦法清楚的告訴我們圖的範圍。剛剛的assignment會把R如何繪製直方圖的細節存至變數x中。請同學輸入:x來顯示繪圖的細節

x

關卡 31

x本質上是一個list,而它名為“breaks”的切點記載了直方圖的切割點。請由$breaks(也就是x$breaks)的值,回答:「樣本數最高的年齡層為:」

(26,28]

關卡 32

這裡幫同學複習:(26,28]的意思是說,age>26且age<=28。x$counts則記載著每個層級的個數。根據x$counts,(26,28]的資料共有45筆。請同學輸入:sum(age>26&age<=28)來驗證。

sum(age > 26 & age <= 28)

關卡 33

hist的結果出發,同學有沒有發現hist(age)本質上,就是把age透過x$breaks切成一段一段,變成順序尺度的數據之後,再用類似剛剛我們繪spon的方式做呈現。

關卡 34

為了接下來的關卡,我們先解釋cut的功能。請同學輸入:?cut

?cut

關卡 35

cut重要的參數有xbreakscut會把breaks當成分割點,將整個實數線分割成若干區段。而資料就會依據區段的位置,被分類為不同的類別。

關卡 36

請同學輸入cut(1:10,2:4)試試看。

cut(1:10, 2:4)

關卡 37

我們可以看到,cut的輸出有三種可能:NA、(2,3]或(3,4]。(2,3]表示原本的數字是大於2小於等於3,(3,4]則表示大於3小於等於4,NA則表示不在前兩種狀況之內。

關卡 38

請同學輸入:plot(cut(age,breaks=x$breaks))

plot(cut(age, breaks = x$breaks))

關卡 39

我們會看到一張圖,以plot(spon)的格式顯示hist(age)的內容。

關卡 40

這種把區點尺度、比例尺度的數據切片後變成順序尺度的技巧,在實務分析上非常常用。如何選取切點,以及切點的個數,在實務上也沒有標準答案,一切都要視其他狀況而定。

關卡 41

舉例來說,我在做政府採購分析時,由於法規會把標案做分類(例如一百萬以下的標案不用公開招標,兩億以上的招標有更多法規需要遵守),此時法規會很自然地告訴我們合理的切點。

關卡 42

histage的例子,則是抓好上下界之後,等距切割。

關卡 43

另外一種做法,是透過quantile函數做切割。這樣每個切片中的資料個數會非常接近,在做進階分析時有好處。

關卡 44

請同學輸入:infert$education,這筆數據本身就是把比例尺度切片後得到順序尺度的結果。

infert$education

關卡 45

在R之中,處理區間尺度與比例尺度的資料,還有其他的方法。請同學輸入:plot(density(age))。

plot(density(age))

關卡 46

這種做法會得到一個平滑的曲線,因為R並沒有先把age做切片,而是用一種叫做kerneldensityestimation的統計方法來估計age這筆數據的probabilitydensityfunction。這個函數在實務上,和長條分配圖一樣是拿來看數據的分佈。

關卡 47

所以我們也可以從這個圖中來觀察:有多少個尖峰?數據在右邊或左邊的尾巴有沒有拉得很長?有沒有明顯的異常值?

關卡 48

請問同學,圖中的曲線,共有幾個尖峰?

2

關卡 49

請問同學,圖中的曲線,是A:左偏(少數特別小的age)、B:右偏(少數特別大的age)或C:不偏(沒有上述兩種狀況)?

C

關卡 50

請問同學,這個曲線有沒有在某個地方有奇怪的尖峰?

No

關卡 51

這個統計方法中,很重要的參數是:bw。我們先請同學來畫出不同bw的結果。請同學先輸入:plot(density(age,bw=0.1))

plot(density(age, bw = 0.1))

關卡 52

bw特別小的狀態,我們會看到一個很醜,而且很接近長條圖的結果。這是因為這筆數據是整數資料,所以當bw小於1很多時,就會看到一個個的尖峰,而這些尖峰的高度則和寬度有關。舉例來說,位在25的尖峰,他的面積必須接近:sum(age==25)/length(age),也就是age為25的樣本數比率。所以越小的bw,尖峰的寬度越小,尖峰就要越高。

關卡 53

請同學輸入:plot(density(age,bw=1))

plot(density(age, bw = 1))

關卡 54

我們可以看到整體的曲線變得平滑許多。但是和plot(density(age))的結果相比,圖形上只剩下兩個峰。同學應該可以就這個例子注意到,bw的參數會影響到結果的拓樸特質,例如峰的個數。

關卡 55

那我們應該如何挑選bw的值呢?這個問題已經有許多統計學家做過相關研究,並且已經有許多學術論文在探討這個問題。而R的density函數的文件中載明了,預設的bw是考量到歷史包袱,所以採用bw.nrd0函數的計算結果。而統計文獻中,比較推薦的方法是bw.SJ函數的計算結果。請同學輸入:plot(density(age,bw="SJ")),使用主流的推薦選bw的方法來繪圖。

plot(density(age, bw = "SJ"))

關卡 56

如果採用nr.SJ方法所選出來的bw參數,得到的結果是只有兩個峰值的。

關卡 57

這裡總結一下上述的討論給同學。我們介紹了一種統計方法來從一筆數值向量中,觀察數據分佈的方法。這裡要強調,我們的目的通常是要觀察數據是否單峰、是否有少數值特別大或特別小,以及不尋常的尖峰。最後我們簡單探討,有些繪圖的參數,是會影響到我們對上述問題的答案。也介紹了R如何實作文獻上的方法,讓我們可以直接站在巨人的肩膀上。

關卡 58

接下來,我們請同學輸入sunspot.year。這是1700年至1988年,每年的太陽黑子數量。

sunspot.year

關卡 59

請問同學,sunspot.year的資料型態是?

class(sunspot.year)

關卡 60

“ts”型態是R專門用來處理時間序列的型態。請同學輸入:plot(sunspot.year)

plot(sunspot.year)

關卡 61

R在繪製時間序列的圖形時,會用線段依序將數據給連接起來。透過這種圖,我們可以判斷資料的值與順序(時間)的相關性強不強。然而,我認為這張圖太密了,資料太多,反而讓我們不容易看出太陽黑子的週期性。請同學輸入:x<-tail(sunspot.year,100),直接取出最後的100筆數據。

x <- tail(sunspot.year, 100)

關卡 62

請同學先檢查x的型態。經過tail函數處理後,x變成數值向量了。

class(x)

關卡 63

請同學還是先輸入:plot(x)

plot(x)

關卡 64

接著再請同學輸入:lines(x)

lines(x)

關卡 65

同學注意到,這次的圖中同時有折線,也有資料點。而且我們也更清楚地看到圖上的資料。

關卡 66

R之中,有非常多修飾圖的函數,例如剛剛同學使用的lines就是其一。當lines收到單一的數值向量後,就會有類似plot(x,type="l")的效果。而lines這類函數不會像plot「另起新圖」,而是在原本的圖上做修飾。因此,在R中稱呼plothist這種會「另起新圖」的函數為「高階繪圖函數」,而lines或是points等做修飾的函數,稱為「低階繪圖函數」。

關卡 67

請同學輸入:lines(x,lty=3,lwd=3,col=2)

lines(x, lty = 3, lwd = 3, col = 2)

關卡 68

這裡的lty參數讓我們可以把線段從直線改成虛線、lwd參數讓我們調整線段的粗細、col讓線段變成紅色。R有大量的低階繪圖函數,讓我們可以對圖形做細膩的調教。然而以現代的眼光來看,R預設的圖形可能不夠精緻了。

關卡 69

R社群的使用者則提供非常多的繪圖套件,想要提供更好的預設視覺效果。例如之後會仔細介紹的ggplot2

關卡 70

接下來我們介紹,要如何輸出R的圖片。

關卡 71

輸出圖片最簡單的方式,就是先用互動式的console畫圖,然後用savePlot將圖用給定的格式儲存到給定的檔名。我們先輸入dst<-tempfile(fileext=".png")來產生一個隨機的暫存檔案名稱。

dst <- tempfile(fileext = ".png")

關卡 72

接著請輸入:savePlot(dst,type="png"),R就會把我們視窗上看到的圖片輸出至檔案。然而在Rstudio環境中savePlot函數不一定能使用。所以我們這邊就不帶著同學練習。

關卡 73

另一種做法是透過設定繪圖引擎(在R中稱為device)。R預設的繪圖引擎是把圖呈現在使用者面前。而我們可以更改繪圖引擎為「把圖以特定格式繪製到檔案中」。舉例來說,png函數就會讓R建立一個繪圖引擎。請同學輸入:png(dst)開啟png的繪圖引擎。

png(dst)

關卡 74

請同學輸入:plot(x)。同學應該可以注意到,右邊的圖像沒有反應。

plot(x)

關卡 75

最後請同學輸入:dev.off(),關閉png的繪圖引擎,讓R把圖片的資料寫入檔案。

try(dev.off(), silent = TRUE)

關卡 76

最後,使用rstudio的同學等等可以輸入:library(rstudioapi)後使用viewer(dst)看看這個圖片。沒有使用rstudio的同學請輸入dst取得儲存路徑後,自行透過作業系統打開該檔案看看繪圖的結果。透過這樣的動作,確認圖片確實寫入這個檔案。

關卡 77

接下來讓我們進入大魔王時間。保險起見,請同學先輸入skip()讓我們重新初始化繪圖環境。

try(dev.off(), silent = TRUE)

關卡 78

hsb是我們事先準備好的資料,裡面記載了學生的基礎資料與學生在各個科目的考試成績。本資料包含11個欄位、200筆記錄(200位學生),第一欄(id)代表該名學生的編號;第二欄(sex)代表學生性別;第三欄(race)代表學生的種族;第四欄(ses)代表學生家庭社經等級;第五欄(schtyp)代表是學校是公立或私立;第六欄(prog)代表學校的類型;第七欄至第十一欄則是各科考試成績。請同學輸入View(hsb)先看看數據。

View(hsb)

關卡 79

請問同學,write欄位應該用哪一個R型態處理?類別請選A,數值請選B

B

關卡 80

請同學用plotdensity的函數組合繪製hsb$write的數據。bw請用“SJ”

plot(density(hsb$write, bw = "SJ"))

關卡 81

請問這張圖一共有幾個峰?

2

關卡 82

請問資料有沒有異常的極大值或極小值?

Yes

關卡 83

請問極小值大概是多少?

-100

關卡 84

其實hsb$write的數據中,是以-99來記載遺失的資料。這在社會科學的數據中很常見。而適當的視覺化可以讓我們在分析之前,就注意到這些資料的異常。

關卡 85

請同學對hsb$sex畫出圓餅圖

pie(table(hsb$sex))

關卡 86

最後我們直接用大量的範例繪圖程式碼,讓同學學習各種圖形繪製的細節調教。

關卡 87

我們看看圓餅圖的範例。結束請按submit()

# 先建立sex物件
sex <- table(hsb$sex)

# 預設圖
pie(sex)

# 加上標題
pie(sex, main = "Sex")

# 改顏色
col <- rainbow(length(sex)) # 利用rainbow函數產生若干種不同的色系
col # 這是RGB的值
pie(sex, main = "Sex", col = col)

# 加上比率
pct <- sex / sum(sex) * 100 # 計算百分比
label <- paste0(names(sex), " ", pct, "%") # 產生說明文字
label # 同學可以比較label, names(sex), pct 和 "%" 的結果
pie(sex, labels = label)

if (FALSE) {
  # 3D pie chart
  library(plotrix) # 同學請自行安裝這個套件
  pie3D(sex, labels = label)
}

# 請回到console輸入`submit()`

關卡 88

我們看看長條圖的範例。結束請按submit()

# 先建立math物件
math <- hsb$math
hist(math)

# 改變title
hist(math, main = "Histogram of math!")

# 改變x軸說明
hist(math, xlab = "Value of math")

# 改變y軸說明
hist(math, ylab = "frequency")

# 改變顏色
hist(math, col = "blue")

# 圖標
legend("topright", "test")

# 改變切割點
hist(math, breaks = 2)
hist(math, breaks = 10)
hist(math, breaks = 20)
# 請回到console輸入`submit()`

關卡 89

我們看看機率分佈圖的範例。結束請按submit()

math <- hsb$math
plot(density(math))
math.sj <- density(math, bw = "SJ")
plot(math.sj)
# 線的粗細
plot(math.sj, lwd = 2) # lwd越大越粗
# 線的型態
plot(math.sj, lty = 2) 
if (FALSE) {
  # 以下指令可以畫出lty的數字與畫圖後的結果
  showLty <- function(ltys, xoff = 0, ...) {
    stopifnot((n <- length(ltys)) >= 1)
    op <- par(mar = rep(.5,4)); on.exit(par(op))
    plot(0:1, 0:1, type = "n", axes = FALSE, ann = FALSE)
    y <- (n:1)/(n+1)
    clty <- as.character(ltys)
    mytext <- function(x, y, txt)
      text(x, y, txt, adj = c(0, -.3), cex = 0.8, ...)
    abline(h = y, lty = ltys, ...); mytext(xoff, y, clty)
    y <- y - 1/(3*(n+1))
    abline(h = y, lty = ltys, lwd = 2, ...)
    mytext(1/8+xoff, y, paste(clty," lwd = 2"))
  }
  showLty(1:6)
}
# 線的顏色
plot(math.sj, col = "red")
# 對線之下的面積著色
polygon(math.sj, col = "red") # 這是一個低階繪圖函數
# 標題
plot(math.sj, main = "math")
# x軸標題
plot(math.sj, xlab = "math")
# 請回到console輸入`submit()`

關卡 90

我們看看折線圖的範例。結束請按submit()

x <- tail(sunspot.year, 50) # 只選出最後50筆資料做圖

# 畫出散布圖
plot(x) 
# 將點連接起來
lines(x) # 低階繪圖函數
# 調色
lines(x, col = "red")
# 加粗
lines(x, lwd = 2)
# 改變線的型態
plot(x) # 重新畫圖
lines(x, lty = 3)
# 標題
plot(x, main = "sunspot")
# 刪除x 軸座標
plot(x, xaxt = "n")
# y 軸座標更改
plot(x, yaxt = "n") # 先刪除y軸座標
axis(2, at = seq(10, 200, 10), labels = seq(10, 200, 10))
# 請回到console輸入`submit()`