R筆記–(5)初聲試啼-簡單的資料分析(迴歸分析)

skydome20

2016/04/08

返回主目錄

【如果喜歡,你也可以免費支持 R 系列筆記!】


本篇目錄

  1. 查看資料集
  2. 敘述統計-繪圖
  3. 資料預處理
  4. 迴歸分析
  5. 變異數分析(anova)
  6. 總結

到現在,你已經學到許多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筆,共有五個欄位:

  1. 花萼長度(Sepal.Length):計算單位是公分。(連續)

  2. 花萼寬度(Sepal.Width):計算單位是公分。(連續)

  3. 花瓣長度(Petal.Length) :計算單位是公分。(連續)

  4. 花瓣寬度(Petal.Width):計算單位是公分。(連續)

  5. 品種(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)是否符合下面三個假設:

  1. 常態性(Normality)

  2. 獨立性(Independence)

  3. 變異數同質性(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~