用不同方法在R中進行泰坦尼克號幸存者預測練習

現有數據維度：
PassengerId
survival 生存 0 = No, 1 = Yes
pclass 票類 社會經濟地位，1 = Upper, 2 = Middle, 3 = Lower
sex 性別
Age 年齡
sibsp 兄弟姐妹/配偶在泰坦尼克號上
parch 父母/孩子在泰坦尼克號上
ticket 票號
fare 客運票價
cabin 艙位數量
embarked 始發港 C = Cherbourg, Q = Queenstown, S = Southampton

讀取數據，可知船上共有891名乘客，其中男性577人，女性314人，乘客生還率為38.38%，有年齡信息的乘客平均年齡為29歲，且數據需進行以下預處理：

1、年齡存在缺失，進行中位數填補

2、姓名這一列中間值也許有用，進行提取

3、cabin存在較多的缺失

4、登錄港口存在較少缺失，使用明顯最多的登錄港口進行填充

5、船票開頭的標記本例中不進行處理分析

> titanic<-read.csv('titanic.csv')

 

pclass應該是分類變量，是否幸存也是分類變量但是后續需要使用暫時不做變化，只把Pclass轉換為因子，確實是

> titanic$Pclass <- factor(titanic$Pclass)

中位數填補年齡缺失，登錄港口缺失填補為S

median(titanic$Age,na.rm = T)
titanic$Age[is.na(titanic$Age)] <- median(titanic$Age,na.rm = T)
titanic$Embarked[is.na(titanic$Embarked)] <- 'S'

提取姓名中間的稱謂

#首先需要自定義一個函數截取出來稱謂
> nametitle <- function(x) {
+     temp <- unlist(strsplit(x, ','))[2]#strsplit結果為list，要提取第二列轉換為unlist
+     temp <- unlist(strsplit(temp, ' '))[2]#姓名稱謂之前存在一個空格，空格之后為我們所需要的稱謂
+ return(temp)}

#將姓名列轉換為字符型向量
titanic$Name <- as.character(titanic$Name)

#測試一下

> nametitle(titanic$Name[1])
[1] "Mr."

#對整列應用提取函數，並賦值給新的一列名為title

> titanic$title <- unlist(lapply(titanic$Name, nametitle))

#未知稱謂標記為未知

> titanic$title[!titanic$title %in% c('Miss.', 'Mr.', 'Mrs.', 'Ms.')] <- 'unknow'

#轉換為因子變量

> titanic$title <- factor(titanic$title)

cabin處理，由結果可知大多數人其實都沒有關聯艙位，或者說可能有但是關聯艙位的數據缺失嚴重

> CabinNum <- function(x) {
+     if (is.na(x)) {
+         return(0)
+     } else {
+         return(length(unlist(strsplit(x, ' '))))
+     }
+ }
> titanic$Cabin <- as.character(titanic$Cabin)
> titanic$CabinNum <- unlist(lapply(titanic$Cabin, CabinNum))

 name、cabin、票號變量基本上可以去除

> titanic$Name <- NULL
> titanic$Cabin <- NULL
> titanic$Ticket <- NULL

 目前的數據已經沒缺失，接下來針對分類變量進行啞變量處理

 

> base <- data.frame(predict(dummyVars(~., data = titanic), titanic))
> describe(base)

 

將數據集進行切割，切分為訓練集和測試集

> trainid <- createDataPartition(base$PassengerId, p = 0.75, list = F)
> train <- base[trainid,]
> test <- base[-trainid,]

方法：邏輯回歸

初次建模

> logistic <- glm(Survived ~ .,
+                 data = train[, -1],
+                 family = 'binomial'(link = 'logit'))
+ 
+ summary(logistic)

 選擇sig較高的值再次進行建模

> logistic <- glm(Survived ~ Pclass.1
+                             + Pclass.2
+                             + Age
+                             + SibSp
+                             + title.Mr.,
+                 data = train[, -1],
+                 family = 'binomial'(link = 'logit'))
+ 
+ summary(logistic)

用之前生成的測試集對模型進行評估

> test$predict <- predict(logistic, test, type = 'response')
+ test$predictClass <- NULL
+ test$predictClass[test$predict >= 0.5] <- 1 #邏輯回歸得到的是一個概率值，而如何去划分需要我們指定，本例中大於等於0.5為正例，即：幸存。0.5也可以划分為0，已經測試過這個例子中不影響預測結
+ test$predictClass[test$predict < 0.5] <- 0
+ table(test$Survived, test$predictClass)

> table(test$Survived)

  0   1 
136  84 

由以上結果可知：測試集內的對象而言，136個0預測正確的有119個，84個1預測正確的有66個，模型效果尚可。

可以再調整一下，我們知道女士更容易獲救，而本次建模中female變量沒有進去（將稱謂中的mr變量換成female），再試試

> logistic <- glm(Survived ~ Pclass.1 + Pclass.2 + Sex.female + SibSp,data = train[,-1],family = 'binomial'(link='logit'))
> summary(logistic)

效果沒有原先的好，在考慮共線性的基礎上，可以自己嘗試疊加不同變量，檢驗模型效果，可同時結合建模目的進行變量挑選。 

方法二：k近鄰

> train <- base[trainid,]
> test <- base[- trainid,]

定義標准化函數並且進行標准化#此里中未處理異常值

> normalize <- function(x) {
+     return((x - min(x)) / (max(x) - min(x))) 
+ }
+ train$Age <- normalize(train$Age)
+ train$CabinNum <- normalize(train$CabinNum)
+ train$SibSp <- normalize(train$SibSp)
+ train$Parch <- normalize(train$Parch)
+ train$Fare <- normalize(train$Fare)
+ test$Age <- normalize(test$Age)
+ test$CabinNum <- normalize(test$CabinNum)
+ test$SibSp <- normalize(test$SibSp)
+ test$Parch <- normalize(test$Parch)
+ test$Fare <- normalize(test$Fare)

> library(class)
> knn(train[,-c(1,2)],test[,-c(1,2)],train$Survived,k = 1)
  [1] 1 0 1 0 1 0 0 1 1 0 0 0 1 0 0 0 0 1 0 0 0 0 1 1 1 0 0 1 1 0 1 0 0 1 0 1 1 0 0 1 1 0 0 1 1 1 1 0 0 1
 [51] 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 1 1 1 0 0 1 1 1 1 0 1 1 1 1 1 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0
[101] 0 0 0 1 0 1 0 0 0 0 1 1 1 0 0 0 0 1 0 0 0 0 0 0 1 0 0 1 1 0 1 0 1 1 1 0 0 1 1 0 1 0 0 0 0 1 0 0 0 1
[151] 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 1 1 1 0 0 0 0 0 1 1 1 0 0 0 0 1 0 1 0 1 0 0 0 1 0 0 0 0 0
[201] 0 0 1 0 0 1 1 1 0 0 0 0 0 1 1 1 1 0 0 1
Levels: 0 1
> table(test$Survived, knn(train[, - c(1, 2)], test[, - c(1, 2)], train$Survived, k = 1))
   
      0   1
  0 114  22
  1  25  59

可以通過調節K的數量來查看模型效果

> table(test$Survived, knn(train[, - c(1, 2)], test[, - c(1, 2)], train$Survived, k = 3))
   
      0   1
  0 122  14
  1  24  60
> table(test$Survived, knn(train[, - c(1, 2)], test[, - c(1, 2)], train$Survived, k = 5))
   
      0   1
  0 124  12
  1  21  63
> table(test$Survived, knn(train[, - c(1, 2)], test[, - c(1, 2)], train$Survived, k = 7))
   
      0   1
  0 126  10
  1  19  65
> table(test$Survived, knn(train[, - c(1, 2)], test[, - c(1, 2)], train$Survived, k = 9))
   
      0   1
  0 127   9
  1  22  62

K=7時，學得最好，不過考慮到運算開銷，可以選擇K=3作為模型，因為效果差異不太大。

K為奇數，因為要投票

方法三：朴素貝葉斯，只能作用在分類變量上。使用啞變量處理之前的數據，選擇其中的分類變量

>install.packages(e1071)
> library(e1071)

> trainid<-createDataPartition(titanic$PassengerId,p = 0.75,list = F)
> train<-titanic[trainid,]
> test<-titanic[-trainid,]
> train$Survived <- factor(train$Survived)
+ test$Survived <- factor(test$Survived)

> fit <- naiveBayes(train[, c('Pclass', 'Sex', 'Embarked', 'title')], train$Survived)
> table(test$Survived, predict(fit, test))
   
      0   1
  0 121  22
  1  23  54
#用laplace修正結果，本例中自變量少，修正后對結果無影響，對於多自變量的模型效果較為明顯

> fit <- naiveBayes(train[, c('Pclass', 'Sex', 'Embarked', 'title')],
+                   train$Survived,
+                   laplace = 5)
+ 
+ table(test$Survived, predict(fit, test))
   
      0   1
  0 121  22
  1  23  54

方法四：決策樹，既能做回歸，也能做分類，這是分類‘class’

trainid <- createDataPartition(titanic$PassengerId,p = 0.75,list = F)
train <- titanic[trainid,]
test <- titanic[- trainid,]

fit <- rpart(Survived ~ .,
             data = train[,-1],
             method = 'class')
> fit
n= 671 

node), split, n, loss, yval, (yprob)
      * denotes terminal node

 1) root 671 265 0 (0.60506706 0.39493294)  
   2) title=Mr. 386  64 0 (0.83419689 0.16580311) *
   3) title=Miss.,Mrs.,unknow 285  84 1 (0.29473684 0.70526316)  
     6) Pclass=3 129  64 0 (0.50387597 0.49612403)  
      12) Fare>=24.80835 28   2 0 (0.92857143 0.07142857) *
      13) Fare< 24.80835 101  39 1 (0.38613861 0.61386139)  
        26) Fare>=13.90835 37  18 0 (0.51351351 0.48648649)  
          52) Fare< 14.8729 8   0 0 (1.00000000 0.00000000) *
          53) Fare>=14.8729 29  11 1 (0.37931034 0.62068966) *
        27) Fare< 13.90835 64  20 1 (0.31250000 0.68750000) *
     7) Pclass=1,2 156  19 1 (0.12179487 0.87820513)  
      14) Sex=male 26  13 0 (0.50000000 0.50000000)  
        28) Age>=17 16   3 0 (0.81250000 0.18750000) *
        29) Age< 17 10   0 1 (0.00000000 1.00000000) *
      15) Sex=female 130   6 1 (0.04615385 0.95384615) *
> fancyRpartPlot(fit)#需要用到四個包:rattle

> table(test$Survived, predict(fit, test, type = 'class'))
   
      0   1
  0 127  16
  1  21  56

 模型評估的另一種方式：

混淆矩陣相關知識：

預測准確率=正確預測的正反例數 / 總數

誤分類率 = 錯誤預測的正反例數 / 總數

真陽性率(召回率)：正確預測到的正例的數量/實際上正例的個數，在實際應用中覆蓋率是很重要的指標，

陽性預測值=正確預測到的正例數 / 預測正例總數

真陰性率 = 正確預測到的負例個數 / 實際負例總數

陰性預測值 = 正確預測到的負例個數 / 預測負例總數