本篇目錄
到現在,你已經學到許多R的技巧了。可是學習歸學習,有時候如果沒有實際去操作體會,其效果僅止於「知道」而已,並沒有達到「運用」的境界。
因此,這裡將會運用之前筆記中提到過的技巧(以及介紹新技巧),進行一次簡單、完整的資料分析,帶大家體會整個資料探勘的流程。
1. 查看資料集
這裡使用的資料,是R內建的鳶尾花(iris)資料(來自於datasets
套件)。
先用str()
和head()
,查看資料裡面的狀態:
require(datasets) # source package
str(iris) # check structure of iris
## 'data.frame': 150 obs. of 5 variables:
## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
head(iris, n=6)
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
不難看出,iris的資料筆數為150筆,共有五個欄位:
花萼長度(Sepal.Length):計算單位是公分。(連續)
花萼寬度(Sepal.Width):計算單位是公分。(連續)
花瓣長度(Petal.Length) :計算單位是公分。(連續)
花瓣寬度(Petal.Width):計算單位是公分。(連續)
品種(Species):可分為Setosa,Versicolor和Virginica。(類別)
此外也可以用summary()
,看各個欄位的基本統計資訊:
summary(iris)
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100
## 1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300
## Median :5.800 Median :3.000 Median :4.350 Median :1.300
## Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199
## 3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800
## Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500
## Species
## setosa :50
## versicolor:50
## virginica :50
##
##
##
2. 敘述統計-繪圖
在進行資料分析之前,先仔細觀察資料,看看能不能從裡面找到一些隱藏資訊。
例如,花萼長度(Sepal.Length)和花萼寬度(Sepal.Width),既然都是花萼,可能會有相關,於是畫圖來看看:
#*** 附上三種繪圖系統的程式碼,以ggplot2輸出 ***#
### Base Plotting System
# plot(x=iris$Sepal.Length, y=iris$Sepal.Width,pch=2)
### Lattice
# require(lattice)
# xyplot(Sepal.Width~Sepal.Length, data=iris)
### ggplot2
require(ggplot2)
ggplot(data=iris) +
geom_point(aes(x=Sepal.Length,
y=Sepal.Width)) +
theme_bw()
嗯…好像看不出什麼東西?
那試試看花瓣長度(Petal.Length)和花瓣寬度(Petal.Width):
#*** 附上三種繪圖系統的程式碼,以ggplot2輸出 ***#
### Base Plotting System
# plot(x=iris$Petal.Length, y=iris$Petal.Width,pch=16)
### Lattice
# require(lattice)
# xyplot(Petal.Width~Petal.Length, data=iris)
### ggplot2
require(ggplot2)
ggplot(data=iris) + # 準備畫布
geom_point(aes(x=Petal.Length, # 散布圖
y=Petal.Width)) +
theme_bw() # 改變主題背景成白色
Bingo,可以觀察出來,花瓣長度和寬度之間,存在著線性關係,而且明顯分成兩群(左下角和右上角),推測可能和種類(Species)有關,左下角的資料可能是屬於同一種類的鳶尾花。
為了確認這一點,我們在上面那張圖標上顏色:
#*** 附上三種繪圖系統的程式碼,以ggplot2輸出 ***#
### Base Plotting System
# plot(x=iris$Petal.Length, y=iris$Petal.Width,pch=16)
# d1 <- iris[iris$Species=="versicolor", ]
# points(x=d1$Petal.Length, y=d1$Petal.Width,pch=16, col="green")
# d2 <- iris[iris$Species=="setosa", ]
# points(x=d2$Petal.Length, y=d2$Petal.Width,pch=16, col="red")
# d3 <- iris[iris$Species=="virginica", ]
# points(x=d3$Petal.Length, y=d3$Petal.Width,pch=16, col="blue")
# legend("topleft", pch=16
# legend=c("setosa","versicolor","virginica"),
# col=c("red", "green", "blue")
# )
### Lattice
# require(lattice)
# xyplot(Petal.Width~Petal.Length,
# data=iris,
# pch=16,
# group=Species,
# auto.key=list(space="top",
# columns=3,
# cex.title=1,
# title="Species Labels",
# pch=16)
# )
### ggplot2
require(ggplot2)
ggplot(data=iris) + # 準備畫布
geom_point(aes(x=Petal.Length, # 散布圖
y=Petal.Width,
color=Species)) + # 把不同品種的資料標上顏色
theme_bw() # 改變主題背景成白色
並且看不同種類的鳶尾花,長度和寬度的盒鬚圖:
#*** 附上三種繪圖系統的程式碼,以ggplot2輸出 ***#
### Base Plotting System
# boxplot(Petal.Length~Species, data=iris, xlab="Species", ylab="Petal.Length")
# boxplot(Petal.Width~Species, data=iris, xlab="Species", ylab="Petal.Length")
### Lattice
# require(lattice)
# bwplot(x = Petal.Length~Petal.Width | Species, data = iris)
### ggplot2
require(ggplot2)
qplot(x=Petal.Length,
y=Petal.Width,
data=iris,
geom="boxplot", # graph type is boxplot
color=Species)
3. 資料預處理
資料探勘的分析過程中,「資料預處理」往往是最花時間的(佔整個流程的70~80%)。
根據不同的資料,預處理手法也會不一樣(改變結構、類別轉啞變數、正規化…),而在預處理之中,最常見的莫過於「遺漏值的處理」!
要用R檢查資料裡是否有遺漏值的存在,需要使用is.na()
的函式:
data <- data.frame(x=c(1,2,3,NA,5),
y=c(4,5,3,NA,NA))
data
## x y
## 1 1 4
## 2 2 5
## 3 3 3
## 4 NA NA
## 5 5 NA
is.na(data) # 遺漏值的地方,標註為TRUE (TRUE/FALSE矩陣的型態)
## x y
## [1,] FALSE FALSE
## [2,] FALSE FALSE
## [3,] FALSE FALSE
## [4,] TRUE TRUE
## [5,] FALSE TRUE
table(is.na(data)) # 資料中總共有多少個遺漏值
##
## FALSE TRUE
## 7 3
有遺漏值的資料會影響分析結果,因此我們會採取一些手段,主要可以分為兩類「移除有遺漏值的資料」、「填補遺漏值」:
# 移除有遺漏值的資料,以下兩種方法都可以 #
data[complete.cases(data), ] # 1.使用 complete.cases()
na.omit(data) # 2.或是使用 na.omit()
# 填補遺漏值(用平均數填值) #
data[is.na(data[,"y"]), "y"] <- mean(data[,"y"], na.rm=T)
data
## x y
## 1 1 4
## 2 2 5
## 3 3 3
## 4 NA 4
## 5 5 4
現在回到iris的資料,檢查看看裡面有沒有遺漏值:
table(is.na(iris))
##
## FALSE
## 750
…看來十分完美,沒有遺漏值,所以讓我們繼續下去吧!
4. 迴歸分析
回歸分析是以一個或一組自變數(解釋變數、預測變項,Xi),來預測一個數值性的因變數(依變數、應變數、被預測變項,Y)。
相信大家都知道,簡單迴歸表示只有一個Y;複回歸則允許多個Y存在。
要在R跑線性回歸的模型,要使用函式lm()
(Linear Model): model = lm(Y ~ X1+X2+…+Xk, data=…)
在這裡,我們以Sepal.Length為依變數(Y),以Sepal.Width、Petal.Length、Petal.Width為自變數(X),進行迴歸分析:
model <- lm(formula= Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width,
data=iris)
summary(model)
##
## Call:
## lm(formula = Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width,
## data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.82816 -0.21989 0.01875 0.19709 0.84570
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.85600 0.25078 7.401 9.85e-12 ***
## Sepal.Width 0.65084 0.06665 9.765 < 2e-16 ***
## Petal.Length 0.70913 0.05672 12.502 < 2e-16 ***
## Petal.Width -0.55648 0.12755 -4.363 2.41e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3145 on 146 degrees of freedom
## Multiple R-squared: 0.8586, Adjusted R-squared: 0.8557
## F-statistic: 295.5 on 3 and 146 DF, p-value: < 2.2e-16
從報表中來看,我們可以獲得許多資訊:
Sepal.Length = 1.85600 + 0.65084xSepal.Width + 0.70913xPetal.Length - 0.55648xPetal.Width
根據p-value,三個自變數(X)對Y都表示顯著。
R-squared: 0.8586 ; Adj R-squared: 0.8557,表示模型預測能力不錯。
Residual standard error: 0.3145
然而,當我們建立出一個線性回歸時,必須要確認其殘差(residual)是否符合下面三個假設:
常態性(Normality)
獨立性(Independence)
變異數同質性(Homogeneity of Variance)
故,首先我們要先從回歸模型中找到殘差的值,可以使用names()
函式,查看回歸模型內具有的資訊:
names(model)
## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "xlevels" "call" "terms" "model"
其中,residuals就是指殘差的值(coefficients代表係數),因此我們可以取出來後進行上面三個假設的檢定:
model$residual
## 1 2 3 4 5
## 0.0845842387 0.2100028184 -0.0492514176 -0.2259940935 -0.0804994772
## 6 7 8 9 10
## 0.0228063193 -0.2946837793 -0.0212452413 -0.2249134657 0.0183576405
## 11 12 13 14 15
## 0.1835036110 -0.2921584372 0.0543545524 -0.2329058599 0.6009920509
## 16 17 18 19 20
## 0.1392141315 0.3064591029 0.1402325047 0.3322417692 -0.1259318390
## 21 22 23 24 25
## 0.2369283669 -0.0051998570 -0.1968466935 0.1689568809 -0.5048980249
## 26 27 28 29 30
## 0.1681764266 0.0191380949 0.1136710428 0.2496679547 -0.2619910053
## 31 32 33 34 35
## -0.0969072894 0.4900512908 -0.3324795188 0.0289982272 0.0740059065
## 36 37 38 39 40
## 0.3216617783 0.5554974346 -0.2361477432 -0.2190839857 0.0787547587
## 41 42 43 44 45
## 0.1111457007 0.3921502918 -0.3492514176 0.0653509110 -0.3539363566
## 46 47 48 49 50
## 0.1656510844 -0.2524933009 -0.2201646135 0.0835036110 0.1147516706
## 51 52 53 54 55
## 0.5074791136 0.1049537714 0.3863847037 0.0339766623 0.3943754392
## 56 57 58 59 60
## -0.4460078969 -0.1463080703 -0.3016594803 0.3179951913 -0.3997967395
## 61 62 63 64 65
## -0.0831510084 -0.0521392090 0.4321155802 -0.1972697386 0.0271271504
## 66 67 68 69 70
## 0.4853024172 -0.5648787967 -0.1642161954 0.5557909307 -0.0365741057
## 71 72 73 74 75
## -0.4408410183 0.3085580827 0.1768869993 -0.2434825547 0.3307347790
## 76 77 78 79 80
## 0.4503861332 0.4969007814 0.2918517557 -0.0997950808 0.2263466961
## 81 82 83 84 85
## -0.0005771938 0.0146877361 0.0889067285 -0.3394585584 -0.7648787967
## 86 87 88 89 90
## -0.3695653944 0.3282110955 0.5503238787 -0.3925225451 -0.0961907695
## 91 92 93 94 95
## -0.5005755351 -0.1914402587 0.0830772485 -0.1365757643 -0.2681845932
## 96 97 98 99 100
## -0.4190840070 -0.2983520251 0.1307347790 0.1016446576 -0.1623551132
## 101 102 103 104 105
## -0.5673452231 -0.3725137603 0.2762260566 -0.4128954378 -0.1972124815
## 106 107 108 109 110
## 0.2798336852 -0.8281636850 0.0907121908 0.1056130341 0.0664904332
## 111 112 113 114 115
## 0.0577159260 0.0856598478 0.2598788402 -0.2157848666 -0.1593561462
## 116 117 118 119 120
## -0.0171656678 -0.2070659578 -0.1561009722 0.5387254932 0.0012249512
## 121 122 123 124 125
## 0.1991815486 -0.4401228184 0.3834396551 0.2136643655 -0.1771986994
## 126 127 128 129 130
## 0.0082006308 0.1194938454 -0.1815867823 -0.0808669238 0.1688979224
## 131 132 133 134 135
## 0.4532705646 0.1453420836 -0.0252186578 -0.1601905403 -0.6402373541
## 136 137 138 139 140
## 0.8456961968 -0.4044244213 -0.3721496737 -0.2106735864 0.3657083202
## 141 142 143 144 145
## 0.1908267264 0.6897444400 -0.3725137603 -0.0426448432 0.0453943647
## 146 147 148 149 150
## 0.4839149600 0.3285668674 0.1169701620 -0.4182462955 -0.5234131742
常態性假設
shapiro.test()
函式可以用來檢驗殘差的常態性:
shapiro.test(model$residual)
##
## Shapiro-Wilk normality test
##
## data: model$residual
## W = 0.99559, p-value = 0.9349
由於虛無假設H0:殘差服從常態分配,因為p-value > 0.05,代表不會拒絕H0。
獨立性假設
要檢驗殘差的獨立性,需要使用套件car
中的durbinWatsonTest()
函式:
require(car)
# 因為這個函式會自動去抓模型中的殘差,故這裡放的是模型,而不是殘差的值
durbinWatsonTest(model)
## lag Autocorrelation D-W Statistic p-value
## 1 -0.03992126 2.060382 0.864
## Alternative hypothesis: rho != 0
由於虛無假設H0:殘差間相互獨立,因為p-value > 0.05,代表不會拒絕H0。
變異數同質性假設
要檢驗殘差的變異數同質性,需要使用套件car
中的ncvTest()
函式:
require(car)
# 因為這個函式會自動去抓模型中的殘差,故這裡放的是模型,而不是殘差的值
ncvTest(model)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 4.448612 Df = 1 p = 0.03492962
由於虛無假設H0:殘差變異數具有同質性,因為p-value < 0.05,代表拒絕H0。(這表示上面的線性模型無法使用)
預測
最後,我們建立模型的目的,是要用來預測!
因此,現在我們手上有一筆新的觀測值,只有Sepal.Width、Petal.Length、Petal.Width的資訊,那我們就可以用建好的迴歸模型,預測出Sepal.Length的值,這時使用predict()
函式:
new.iris <- data.frame(Sepal.Width=3.456, Petal.Length=1.535, Petal.Width=0.341)
new.iris
## Sepal.Width Petal.Length Petal.Width
## 1 3.456 1.535 0.341
predict(model, new.iris)
## 1
## 5.004048
5. 變異數分析(anova)
經過視覺化的步驟,發現三個品種鳶尾花的Petal.Width或Petal.Length(平均數)有所差異。
若要用統計上的檢定,要進一步地確認,就可以使用變異數分析(anova)。
假設檢定的對應H0和H1分別如下:
H0: μ(Setosa) = μ(Versicolor) = μ(Virginica)
H1: 至少有一種平均數和其他品種不相等
要用o ne-way-anova ,R的函式是anova()
,並且事先要跑線性迴歸模型:
a.lm <- lm(Petal.Width~Species, data=iris)
anova(a.lm)
## Analysis of Variance Table
##
## Response: Petal.Width
## Df Sum Sq Mean Sq F value Pr(>F)
## Species 2 80.413 40.207 960.01 < 2.2e-16 ***
## Residuals 147 6.157 0.042
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
b.lm <- lm(Petal.Length~Species, data=iris)
anova(b.lm)
## Analysis of Variance Table
##
## Response: Petal.Length
## Df Sum Sq Mean Sq F value Pr(>F)
## Species 2 437.10 218.551 1180.2 < 2.2e-16 ***
## Residuals 147 27.22 0.185
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
兩者的 p-value 都遠小於0.05,表示不同品種間確實有顯著差異。
總結
完成這篇筆記後,除了複習之前的技巧之外,還學到了新的技巧:遺漏值處理,迴歸分析,變異數分析。
事實上,要學會R的各種技巧並不難!難的是當我們陸續學到許多技巧後,要如何把這些技巧靈活運用在各式各樣的資料上。同時,你需要了解R再強大,充其量不過只是一個工具而已,若沒有紮實基礎與清楚的思維,也只是在舞刀耍棍罷了,實際上是派不上用場的。
之後,會繼續介紹各種不同的模型(決策樹、類神經網路…),在R上怎麼操作,並且根據不同的資料,導入不同(資料)處理手法。
It’s still a long way to go~