完整的R語言預測建模實例-從數據清理到建模預測

本文使用Kaggle上的一個公開數據集，從數據導入，清理整理一直介紹到最后數據多個算法建模，交叉驗證以及多個預測模型的比較全過程，注重在實際數據建模過程中的實際問題和挑戰，主要包括以下五個方面的挑戰：

    缺失值的挑戰
    異常值的挑戰
    不均衡分布的挑戰
    (多重)共線性的挑戰
   預測因子的量綱差異
以上的幾個主要挑戰，對于熟悉機器學習的人來說，應該都是比較清楚的，這個案例中會涉及到五個挑戰中的缺失值，量綱和共線性問題的挑戰。
案例數據說明
本案例中的數據可以在下面的網址中下載：
https://www.kaggle.com/primaryobjects/voicegender/downloads/voicegender.zip
下載到本地后解壓縮會生成voice.csv文件
下面首先大概了解一下我們要用來建模的數據

數據共包含21個變量，最后一個變量label是需要我們進行預測的變量，即性別是男或者女
前面20個變量都是我們的預測因子，每一個都是用來描述聲音的量化屬性。
下面我們開始我們的具體過程
步驟1：基本準備工作

步驟1主要包含以下三項工作：

    設定工作目錄
    載入需要使用的包
    準備好并行計算

### the first step: set your working directory
setwd("C:/Users/chn-fzj/Desktop/R Projects/Kaggle-Gender by Voice")
### R中的文件路徑應把Windows系統默認的"\"替換為"/"

### load packages to be used, if not installed, please use ##install.packages("yourPackage")
require(readr)
require(ggplot2)
require(dplyr)
require(tidyr)
require(caret)
require(corrplot)
require(Hmisc)
require(parallel)
require(doParallel)
require(ggthemes)

# parallel processing set up
n_Cores <- detectCores()##檢測你的電腦的CPU核數
cluster_Set <- makeCluster(n_Cores)##進行集群
registerDoParallel(cluster_Set)

步驟2：數據的導入和理解

數據下載解壓縮后就是一份名為‘voice.csv’ 的文件，我們將csv文件存到我們設定的工作目錄之中，就可以導入數據了。

### read in original dataset
voice_Original <- read_csv("voice.csv",col_names=TRUE)
describe(voice_Original)
Hmisc包中的describe 函數是我個人最喜歡的對數據集進行概述，整體上了解數據集的最好的一個函數，運行結果如下：

voice_Original

 21  Variables      3168  Observations
-------------------------------------------------------------------
meanfreq
      n missing  unique    Info    Mean     .05     .10     .25
   3168       0    3166       1  0.1809  0.1260  0.1411  0.1637
    .50     .75     .90     .95
 0.1848  0.1991  0.2177  0.2291

lowest : 0.03936 0.04825 0.05965 0.05978 0.06218
highest: 0.24353 0.24436 0.24704 0.24964 0.25112
-------------------------------------------------------------------
sd
      n missing  unique    Info    Mean     .05     .10     .25
   3168       0    3166       1 0.05713 0.03162 0.03396 0.04195
    .50     .75     .90     .95
0.05916 0.06702 0.07966 0.08549

lowest : 0.01836 0.02178 0.02400 0.02427 0.02456
highest: 0.11126 0.11126 0.11265 0.11451 0.11527
-------------------------------------------------------------------
median
      n missing  unique    Info    Mean     .05     .10     .25
   3168       0    3077       1  0.1856  0.1164  0.1340  0.1696
    .50     .75     .90     .95
 0.1900  0.2106  0.2274  0.2358

lowest : 0.01097 0.01359 0.01579 0.02699 0.02936
highest: 0.25663 0.25698 0.25742 0.26054 0.26122
-------------------------------------------------------------------
Q25
      n missing  unique    Info    Mean     .05     .10     .25
   3168       0    3103       1  0.1405 0.04358 0.07509 0.11109
    .50     .75     .90     .95
0.14029 0.17594 0.20063 0.21524

lowest : 0.0002288 0.0002355 0.0002395 0.0002502 0.0002669
highest: 0.2394595 0.2405416 0.2407352 0.2421235 0.2473469
-------------------------------------------------------------------
Q75
      n missing  unique    Info    Mean     .05     .10     .25
   3168       0    3034       1  0.2248  0.1874  0.1963  0.2087
    .50     .75     .90     .95
 0.2257  0.2437  0.2536  0.2577

lowest : 0.04295 0.05827 0.07596 0.09019 0.09267
highest: 0.26879 0.26892 0.26894 0.26985 0.27347
-------------------------------------------------------------------
IQR
      n missing  unique    Info    Mean     .05     .10     .25
   3168       0    3073       1 0.08431 0.02549 0.02931 0.04256
    .50     .75     .90     .95
0.09428 0.11418 0.13284 0.15632

lowest : 0.01456 0.01492 0.01511 0.01549 0.01659
highest: 0.24530 0.24597 0.24819 0.24877 0.25223
-------------------------------------------------------------------
skew
      n missing  unique    Info    Mean     .05     .10     .25
   3168       0    3166       1    3.14   1.123   1.299   1.650
    .50     .75     .90     .95
  2.197   2.932   3.916   6.918

lowest :  0.1417  0.2850  0.3260  0.5296  0.5487
highest: 32.3507 33.1673 33.5663 34.5375 34.7255
-------------------------------------------------------------------
kurt
      n missing  unique    Info    Mean     .05     .10     .25
   3168       0    3166       1   36.57   3.755   4.293   5.670
    .50     .75     .90     .95
  8.318  13.649  27.294  75.169

lowest :    2.068    2.210    2.269    2.293    2.463
highest: 1128.535 1193.434 1202.685 1271.354 1309.613
-------------------------------------------------------------------
sp.ent
      n missing  unique    Info    Mean     .05     .10     .25
   3168       0    3166       1  0.8951  0.8168  0.8322  0.8618
    .50     .75     .90     .95
 0.9018  0.9287  0.9513  0.9630

lowest : 0.7387 0.7476 0.7477 0.7485 0.7487
highest: 0.9764 0.9765 0.9765 0.9785 0.9820
-------------------------------------------------------------------
sfm
      n missing  unique    Info    Mean     .05     .10     .25
   3168       0    3166       1  0.4082  0.1584  0.1883  0.2580
    .50     .75     .90     .95
 0.3963  0.5337  0.6713  0.7328

lowest : 0.03688 0.08024 0.08096 0.08220 0.08266
highest: 0.82259 0.82267 0.82610 0.83135 0.84294
-------------------------------------------------------------------
mode
      n missing  unique    Info    Mean     .05     .10     .25
   3168       0    2825       1  0.1653 0.00000 0.01629 0.11802
    .50     .75     .90     .95
0.18660 0.22110 0.24901 0.26081

lowest : 0.0000000 0.0007279 0.0007749 0.0008008 0.0008427
highest: 0.2791181 0.2795230 0.2795852 0.2797034 0.2800000
-------------------------------------------------------------------
centroid
      n missing  unique    Info    Mean     .05     .10     .25
   3168       0    3166       1  0.1809  0.1260  0.1411  0.1637
    .50     .75     .90     .95
 0.1848  0.1991  0.2177  0.2291

lowest : 0.03936 0.04825 0.05965 0.05978 0.06218
highest: 0.24353 0.24436 0.24704 0.24964 0.25112
-------------------------------------------------------------------
meanfun
      n missing  unique    Info    Mean     .05     .10     .25
   3168       0    3166       1  0.1428 0.09363 0.10160 0.11700
    .50     .75     .90     .95
0.14052 0.16958 0.18519 0.19343

lowest : 0.05557 0.05705 0.06097 0.06254 0.06348
highest: 0.22342 0.22576 0.22915 0.23114 0.23764
-------------------------------------------------------------------
minfun
      n missing  unique    Info    Mean     .05     .10     .25
   3168       0     913       1  0.0368 0.01579 0.01613 0.01822
    .50     .75     .90     .95
0.04611 0.04790 0.05054 0.05644

lowest : 0.009775 0.009785 0.009901 0.009911 0.010163
highest: 0.168421 0.178571 0.185185 0.200000 0.204082
-------------------------------------------------------------------
maxfun
      n missing  unique    Info    Mean     .05     .10     .25
   3168       0     123    0.99  0.2588  0.1925  0.2192  0.2540
    .50     .75     .90     .95
 0.2712  0.2775  0.2791  0.2791

lowest : 0.1031 0.1053 0.1087 0.1111 0.1124
highest: 0.2774 0.2775 0.2778 0.2791 0.2791
-------------------------------------------------------------------
meandom
      n missing  unique    Info    Mean     .05     .10     .25
   3168       0    2999       1  0.8292  0.1045  0.1888  0.4198
    .50     .75     .90     .95
 0.7658  1.1772  1.5602  1.8004

lowest : 0.007812 0.007979 0.007990 0.008185 0.008247
highest: 2.544271 2.591580 2.676989 2.805246 2.957682
-------------------------------------------------------------------
mindom
       n  missing   unique     Info     Mean      .05      .10
    3168        0       77     0.92  0.05265 0.007812 0.007812
     .25      .50      .75      .90      .95
0.007812 0.023438 0.070312 0.164062 0.187500

lowest : 0.004883 0.007812 0.014648 0.015625 0.019531
highest: 0.343750 0.351562 0.400391 0.449219 0.458984
-------------------------------------------------------------------
maxdom
      n missing  unique    Info    Mean     .05     .10     .25
   3168       0    1054       1   5.047  0.3125  0.6094  2.0703
    .50     .75     .90     .95
 4.9922  7.0078  9.4219 10.6406

lowest :  0.007812  0.015625  0.023438  0.054688  0.070312
highest: 21.515625 21.562500 21.796875 21.843750 21.867188
-------------------------------------------------------------------
dfrange
      n missing  unique    Info    Mean     .05     .10     .25
   3168       0    1091       1   4.995  0.2656  0.5607  2.0449
    .50     .75     .90     .95
 4.9453  6.9922  9.3750 10.6090

lowest :  0.000000  0.007812  0.015625  0.019531  0.024414
highest: 21.492188 21.539062 21.773438 21.820312 21.843750
-------------------------------------------------------------------
modindx
      n missing  unique    Info    Mean     .05     .10     .25
   3168       0    3079       1  0.1738 0.05775 0.07365 0.09977
    .50     .75     .90     .95
0.13936 0.20918 0.32436 0.40552

lowest : 0.00000 0.01988 0.02165 0.02194 0.02217
highest: 0.84448 0.85470 0.85776 0.87950 0.93237
-------------------------------------------------------------------
label
      n missing  unique
   3168       0       2

female (1584, 50%), male (1584, 50%)
-------------------------------------------------------------------
通過這個函數，我們現在可以對數據集中的每一個變量都有一個整體性把握。
我們可以看出我們共有21個變量，共計3168個觀測值。

由于本數據集數據完整，沒有缺失值，因而我們實際上并沒有缺失值的挑戰，但是為了跟實際的數據挖掘過程相匹配，我們會人為將一些數據設置為缺失值，并對這些缺失值進行插補，大家也可以實際看一下我們應用的插補法的效果：

###missing values
## set 30 numbers in the first column into NA
set.seed(1001)
random_Number <- sample(1:3168,30)
voice_Original1 <- voice_Original
voice_Original[random_Number,1] <- NA
describe(voice_Original)

meanfreq
      n missing  unique    Info    Mean     .05     .10     .25
   3138      30    3136       1  0.1808  0.1257  0.1411  0.1635
    .50     .75     .90     .95
 0.1848  0.1991  0.2176  0.2291

lowest : 0.03936 0.04825 0.05965 0.05978 0.06218
highest: 0.24353 0.24436 0.24704 0.24964 0.25112

這時候我們能看見，第一個變量meanfreq 中有了30個缺失值，現在我們需要對他們進行插補，我們會用到caret 包中的preProcess 函數

### impute missing data
original_Impute <- preProcess(voice_Original,method="bagImpute")
voice_Original <- predict(original_Impute,voice_Original)

現在我們來看一下我們插補法的結果，我們的方法就是將我們設為缺失值的原始值和我們插補后的值結合到一個數據框中，大家可以進行直接比較：

### compare results of imputation
compare_Imputation <- data.frame(
  voice_Original1[random_Number,1],
  voice_Original[random_Number,1]
)
compare_Imputation

對比結果如下：

    meanfreq meanfreq.1
1  0.2122875  0.2117257
2  0.1826562  0.1814900
3  0.2009399  0.1954627
4  0.1838745  0.1814900
5  0.1906527  0.1954627
6  0.2319645  0.2313031
7  0.1736314  0.1814900
8  0.2243824  0.2313031
9  0.1957448  0.1954627
10 0.2159557  0.2117257
11 0.2047696  0.2084277
12 0.1831099  0.1814900
13 0.1873643  0.1814900
14 0.2077344  0.2117257
15 0.1648246  0.1651041
16 0.1885224  0.1898701
17 0.1342805  0.1272604
18 0.1933665  0.1954627
19 0.1888149  0.1940667
20 0.2180404  0.2117257
21 0.1980392  0.1954627
22 0.1898704  0.1954627
23 0.1761953  0.1814900
24 0.2356528  0.2313031
25 0.1785359  0.1814900
26 0.1856824  0.1814900
27 0.1808664  0.1814900
28 0.1784912  0.1814900
29 0.1990789  0.1954627
30 0.1714903  0.1651041

可以看出，我們的插補出來的值和原始值之間的差異是比較小的，可以幫助我們進行下一步的建模工作。

另外一點，我們在實際工作中，我們用到的預測因子中，往往包含數值型和類別型的數據，但是我們數據中全部都是數值型的，所以我們要增加難度，將其中的一個因子轉換為類別型數據，具體操作如下：

### add a categorcial variable
voice_Original <- voice_Original%>%
  mutate(sp.ent=
           ifelse(sp.ent>0.9,"High","Low"))
除了使用describe 函數掌握數據的基本狀況外，一個必不可少的數據探索步驟，就是使用圖形進行探索，我們這里只使用一個例子，幫助大家了解：

### visual exploration of the dataset
voice_Original%>%
  ggplot(aes(x=meanfreq,y=dfrange))+
  geom_point(aes(color=label))+
  theme_wsj()


圖形結果如下：

但是我們更關注的是，預測因子之間是不是存在高度的相關性，因為預測因子間的香瓜性對于一些模型，是有不利的影響的。
對于研究預測因子間的相關性，corrplot 包中的corrplot函數提供了很直觀的圖形方法：

###find correlations between factors
factor_Corr <- cor(voice_Original[,-c(9,21)])
corrplot(factor_Corr,method="number")

這個相關性矩陣圖可以直觀地幫助我們發現因子間的強相關性。
步驟3：數據分配與建模

在實際建模過程中，我們不會將所有的數據全部用來進行訓練模型，因為相比較模型數據集在訓練中的表現，我們更關注模型在訓練集，也就是我們的模型沒有遇到的數據中的預測表現。
因此，我們將我們的數據集的70%的數據用來訓練模型，剩余的30%用來檢驗模型預測的結果。

### separate dataset into training and testing sets
sample_Index <- createDataPartition(voice_Original$label,p=0.7,list=FALSE)
voice_Train <- voice_Original[sample_Index,]
voice_Test <- voice_Original[-sample_Index,]

但是我們還沒有解決之前我們發現的一些問題，數據的量綱實際上是不一樣的，另外某些因子間存在高度的相關性，這對我們的建模是不利的，因此我們需要進行一些預處理，我們又需要用到preProcess 函數：

### preprocess factors for further modeling
pp <- preProcess(voice_Train,method=c("scale","center","pca"))
voice_Train <- predict(pp,voice_Train)
voice_Test <- predict(pp,voice_Test)

我們首先將數值型因子進行了標準化，確保所有的因子在一個量綱上，接著對已經標準化的數據進行主成分分析，消除因子中的高相關性。如果我們看一下我們的現在經過處理的數據，就可以看到：

voice_Train

 12  Variables      2218  Observations
-----------------------------------------------------------
sp.ent
      n missing  unique
   2218       0       2

High (1144, 52%), Low (1074, 48%)
-----------------------------------------------------------
label
      n missing  unique
   2218       0       2

female (1109, 50%), male (1109, 50%)
-----------------------------------------------------------
PC1
        n   missing    unique      Info      Mean
     2218         0      2216         1 2.084e-17
      .05       .10       .25       .50       .75
  -5.2623   -3.8212   -2.0470    0.3775    2.0260
      .90       .95
   3.6648    4.5992

lowest : -9.885 -9.138 -8.560 -8.476 -8.412
highest:  6.377  6.381  6.391  6.755  6.934
-----------------------------------------------------------
PC2
         n    missing     unique       Info       Mean
      2218          0       2216          1 -4.945e-16
       .05        .10        .25        .50        .75
   -2.7216    -2.0700    -0.8694     0.2569     0.9934
       .90        .95
    1.5576     2.0555

lowest : -5.528 -5.315 -5.132 -5.103 -5.019
highest:  4.493  4.509  4.598  4.732  4.931
-----------------------------------------------------------
PC3
        n   missing    unique      Info      Mean
     2218         0      2216         1 1.579e-16
      .05       .10       .25       .50       .75
  -1.6818   -1.3640   -0.7880   -0.2214    0.5731
      .90       .95
   1.1723    1.6309

lowest : -2.809 -2.536 -2.462 -2.443 -2.407
highest:  8.055  8.299  8.410  8.805  9.229
-----------------------------------------------------------
PC4
         n    missing     unique       Info       Mean
      2218          0       2216          1 -3.583e-16
       .05        .10        .25        .50        .75
  -1.98986   -1.60536   -0.75468    0.09347    0.86320
       .90        .95
   1.49494    1.83657

lowest : -7.887 -6.616 -5.735 -5.568 -4.596
highest:  2.888  2.921  3.046  3.123  3.311
-----------------------------------------------------------
PC5
         n    missing     unique       Info       Mean
      2218          0       2216          1 -1.127e-16
       .05        .10        .25        .50        .75
   -1.8479    -1.2788    -0.5783     0.0941     0.6290
       .90        .95
    1.1909     1.5739

lowest : -4.595 -3.900 -3.887 -3.787 -3.760
highest:  3.160  3.313  3.548  3.722  3.822
-----------------------------------------------------------
PC6
        n   missing    unique      Info      Mean
     2218         0      2216         1 6.421e-18
      .05       .10       .25       .50       .75
 -1.56253  -1.03095  -0.39648   0.03999   0.53475
      .90       .95
  1.10113   1.38224

lowest : -6.971 -6.530 -5.521 -5.510 -5.320
highest:  1.943  1.948  2.005  2.053  2.066
-----------------------------------------------------------
PC7
         n    missing     unique       Info       Mean
      2218          0       2216          1 -2.789e-16
       .05        .10        .25        .50        .75
   -1.0995    -0.8375    -0.4970    -0.1234     0.4493
       .90        .95
    1.1055     1.4462

lowest : -3.370 -3.132 -2.977 -2.813 -2.664
highest:  2.951  3.136  3.863  3.937  4.128
-----------------------------------------------------------
PC8
         n    missing     unique       Info       Mean
      2218          0       2216          1 -7.291e-17
       .05        .10        .25        .50        .75
  -1.18707   -0.96343   -0.51065   -0.02345    0.46939
       .90        .95
   0.96676    1.28817

lowest : -2.644 -2.611 -2.477 -2.328 -2.261
highest:  2.926  2.940  2.967  2.971  3.456
-----------------------------------------------------------
PC9
        n   missing    unique      Info      Mean
     2218         0      2216         1 4.008e-16
      .05       .10       .25       .50       .75
 -1.06437  -0.84861  -0.47079  -0.04825   0.42092
      .90       .95
  0.96161   1.25187

lowest : -2.267 -2.263 -2.095 -2.066 -1.898
highest:  2.217  2.244  2.266  2.414  2.460
-----------------------------------------------------------
PC10
        n   missing    unique      Info      Mean
     2218         0      2216         1 2.387e-16
      .05       .10       .25       .50       .75
 -0.93065  -0.71784  -0.40541  -0.07025   0.37068
      .90       .95
  0.82534   1.12412

lowest : -2.160 -1.810 -1.754 -1.744 -1.661
highest:  2.164  2.292  2.349  2.385  2.654
-----------------------------------------------------------

原來的所有數值型因子已經被PC1-PC10取代了。

現在，我們進行一些通用的設置，為不同的模型進行交叉驗證比較做好準備。

### define formula
model_Formula <- label~PC1+PC2+PC3+PC4+PC5+PC6+PC7+PC8+PC9+PC10+sp.ent


###set cross-validation parameters

modelControl <- trainControl(method="repeatedcv",number=5,
                             repeats=5,allowParallel=TRUE)

下面我們開始建立我們的第一個模型：邏輯回歸模型：

### model 1: logistic regression
glm_Model <- train(model_Formula,
                   data=voice_Train,
                   method="glm",
                   trControl=modelControl)

將模型應用到測試集上，并將結果與真實值進行比較：

voice_Test1 <- voice_Test[,-2]
voice_Test1$glmPrediction <- predict(glm_Model,voice_Test1)
table(voice_Test$label,voice_Test1$glmPrediction)

我們得到的預測結果如下：

         female male
  female    459   16
  male        7  468

我們的邏輯回歸你模型將7個女性錯判成了男性，16個男性錯判成了女性，應該說結果還是不錯的。
下面我們再來看看下一個模型：線性判別分析(LDA):

### model 2:linear discrimant analysis
lda_Model <- train(model_Formula,
                   data=voice_Train,
                   method="lda",
                   trControl=modelControl)

voice_Test1$ldaPrediction <- predict(lda_Model,voice_Test1)
table(voice_Test$label,voice_Test1$ldaPrediction)

         female male
  female    454   21
  male        6  469
目前lda方法的預測結果略差于邏輯回歸；
第三個模型：隨機森林

### model 3: random forrest
rf_Model <- train(model_Formula,
                  data=voice_Train,
                  method="rf",
                  trControl=modelControl,
                  ntrees=500)
voice_Test1$rfPrediction <- predict(rf_Model,voice_Test1)
table(voice_Test$label,voice_Test1$rfPrediction)

         female male
  female    457   18
  male        6  469

可以看到隨機森林的結果介于上面兩個模型之間。但是模型的結果是存在一定的偶然性的，即因為都使用了交叉驗證，每個模型都存在抽樣的問題，因此結果之間存在一定的偶然性，所以我們需要對模型進行統計意義上的比較。
但是在此之前，我想提一下并行計算的問題，我們在開始建模之前就使用parallel 和doParallel 兩個包設置了并行計算的參數，在modelControl中將allowParallel的值設為了TRUE,就可以幫助我們進行交叉驗證時進行并行計算，下面這張圖可以幫助我們看到差異：

因為原生的R只支持單進程，通過我們的設置，可以將四個核都使用起來，可以大為減少我們的計算時間。

我們最后的一個步驟就是要將三個模型進行比較，確定我們最優的一個模型：

### which model is the best?
model_Comparison <-
  resamples(list(
    LogisticRegression=glm_Model,
    LinearDiscrimant=lda_Model,
    RandomForest=rf_Model
  ))

summary(model_Comparison)
bwplot(model_Comparison,layout=c(2,1))
下面是我們比較的結果：

Call:
summary.resamples(object = model_Comparison)

Models: LogisticRegression, LinearDiscrimant, RandomForest
Number of resamples: 25

Accuracy
                     Min. 1st Qu. Median   Mean 3rd Qu.
LogisticRegression 0.9572  0.9640 0.9685 0.9699  0.9752
LinearDiscrimant   0.9550  0.9640 0.9662 0.9677  0.9729
RandomForest       0.9505  0.9595 0.9640 0.9641  0.9685
                     Max. NA's
LogisticRegression 0.9819    0
LinearDiscrimant   0.9842    0
RandomForest       0.9774    0

Kappa
                     Min. 1st Qu. Median   Mean 3rd Qu.
LogisticRegression 0.9144  0.9279 0.9369 0.9398  0.9505
LinearDiscrimant   0.9099  0.9279 0.9324 0.9354  0.9457
RandomForest       0.9009  0.9189 0.9279 0.9282  0.9369
                     Max. NA's
LogisticRegression 0.9639    0
LinearDiscrimant   0.9685    0
RandomForest       0.9549    0

結果從準確率和Kappa值兩個方面對數據進行了比較，可以幫助我們了解模型的實際表現，當然我們也可以通過圖形展現預測結果：

根據結果，我們可以看到，其實邏輯回歸的結果還是比較好的。
所以我們可以將邏輯回歸的結果作為我們最終使用的模型。