簡介：泰坦尼克號（RMS Titanic，又稱鐵達尼號）是一艘奧林匹克級郵輪，於1912年4月處女航時撞上冰山後沉沒。泰坦尼克號由位於愛爾蘭島貝爾法斯特的哈蘭德與沃爾夫造船廠興建，是當時最大的客運輪船。在它的處女航中，泰坦尼克號從英國南安普敦出發，途經法國瑟堡-奧克特維爾以及愛爾蘭昆士敦，計劃中的目的地為美國紐約。1912年4月14日，船上時間夜裡11時40分，泰坦尼克號撞上冰山；4月15日凌晨2時20分，船體斷裂成兩截後沉入大西洋。泰坦尼克號海難為和平時期死傷人數最慘重的海難之一。船上1500多人喪生。

kaggle將此數據作為預測生存率數據集。以下為整理建模預測過程

首先下載數據-------泰坦尼克號數據集

#載入包及導入數據nnlibrary(ggplot2) # 可視化nlibrary(ggthemes) # 可視化nlibrary(scales) # 可視化nlibrary(dplyr) # 數據處理nlibrary(mice) # 缺失值處理nlibrary(randomForest) # 建模nlibrary(Amelia)#缺失值nlibrary(ipred)#建模n##載入數據ntrain <- read.csv(C:/Users/ximen/Documents/MYR/P5/train.csv,n stringsAsFactors = F,header=T,na.strings=c(""))ntest <- read.csv(C:/Users/ximen/Documents/MYR/P5/test.csv,n stringsAsFactors = F,header=T,na.strings=c(""))nfull <- bind_rows(train, test) # 將訓練和測試集合併n

#查看數據集n> str(full)ndata.frame:t1309 obs. of 20 variables:n $ PassengerId: Factor w/ 1309 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...n $ Survived : int 0 1 1 1 0 0 0 0 1 1 ...n $ Pclass : Factor w/ 3 levels "1","2","3": 3 1 3 1 3 3 1 3 3 2 ...n $ Name : chr "Braund, Mr. Owen Harris" "Cumings, Mrs. John Bradley (Florence Briggs Thayer)" "Heikkinen, Miss. Laina" "Futrelle, Mrs. Jacques Heath (Lily May Peel)" ...n $ Sex : Factor w/ 2 levels "female","male": 2 1 1 1 2 2 2 2 1 1 ...n $ Age : num 22 38 26 35 35 27 54 2 27 14 ...n $ SibSp : int 1 1 0 1 0 0 0 3 0 1 ...n $ Parch : int 0 0 0 0 0 0 0 1 2 0 ...n $ Ticket : chr "A/5 21171" "PC 17599" "STON/O2. 3101282" "113803" ...n $ Fare : num 7.25 71.28 7.92 53.1 8.05 ...n $ Cabin : chr NA "C85" NA "C123" ...n $ Embarked : Factor w/ 3 levels "C","Q","S": 3 1 3 3 3 2 3 3 3 1 ...n

nnSurvived ：n取值0和1（0表示死亡，1表示獲救）

nPclass ：乘客的船倉等級（1，2，3）

nName ：乘客姓名

nSex ：乘客性別

nAge ：乘客年齡

nSibSp ：船上配偶或兄妹的人數

nParch ：船上父母或孩子的人數

nTicket ：票號

nFare ：票價

nCabin ：乘客船艙號

nEmbarked ：出發港口

#從變數-名稱列里提取名稱抬頭nfull$Title <- gsub((.*, )|(..*), , full$Name)n# 展示名稱及性別列聯表ntable(full$Sex, full$Title)n# 將比較少出現的名稱抬頭合併為一個組nrare_title <- c(Dona, Lady, the Countess,Capt, Col, Don, n Dr, Major, Rev, Sir, Jonkheer)n# 修改下部分抬頭名稱nfull$Title[full$Title == Mlle] <- Miss nfull$Title[full$Title == Ms] <- Missnfull$Title[full$Title == Mme] <- Mrs nfull$Title[full$Title %in% rare_title] <- Rare Titlenn# 再觀察下名稱及性別列聯表nn> table(full$Sex, full$Title)n n Master Miss Mr Mrs Rare Titlen female 0 264 0 198 4n male 61 0 757 0 25n

#提取乘客名字nfull$Surname <- sapply(full$Name, n function(x) strsplit(x, split = [,.])[[1]][1])n# 創建一個家庭成員大小變數（數量包含乘客本身）nfull$Fsize <- full$SibSp + full$Parch + 1n# 創建家庭變數，姓氏_家庭人數nfull$Family <- paste(full$Surname, full$Fsize, sep=_)n# 觀察家庭大小及成員生存狀態是否存在相關nggplot(full[1:891,], aes(x = Fsize, fill = factor(Survived))) +n geom_bar(stat=count, position=dodge) +n scale_x_continuous(breaks=c(1:11)) +n labs(x = Family Size) +ntheme_few()n n

# 將家庭大小離散化nfull$FsizeD[full$Fsize == 1] <- singletonnfull$FsizeD[full$Fsize < 5 & full$Fsize > 1] <- smallnfull$FsizeD[full$Fsize > 4] <- largenn# 展示家庭大小及生存關係的馬賽克圖nmosaicplot(table(full$FsizeD, full$Survived), n main=Family Size by Survival, shade=TRUE)n

缺失值處理：

# 查看數據集缺失數據nmissmap(full, main = "Missing values vs observed")n#具體缺失參數，空值nsapply(full,function(x) sum(is.na(x)))nnfull$Deck<-factor(sapply(full$Cabin, function(x) strsplit(x, NULL)[[1]][1]))n

變數Cabin和age有較多的缺失值，survived的缺失值是由於合併了測試集而產生的。用sapply函數顯示full數據集各個變數缺失值情況及變數的因子長度

數據集的數量有限，應盡量保持數據集的完整性。從上可看出各變數缺失值個數，Embarked有2個缺失值，年齡有263個缺失值，Fare有1個缺失值，Cabin變數暫時忽略不使用，survived缺失是因測試數據集引起不必理會。接下來將對這3個變數進行缺失值填補。

1,處理embarked缺失

#處理embarked變數缺失值,首先確定缺失值觀測值所在行數nwhich(full$Embarked %in% NA)n#輸出缺失值信息nfull[full$PassengerId==62,]nfull[full$PassengerId==830,]n

從上面輸出的結果可以發現，這兩名乘客都屬於第一類等級客戶，支付的票價都是80美金。票價與登船港口根據我們的生活經驗可知是具有相關性的。當我們乘坐高鐵時，也是如此，在不同的車站出發，以不同的價格，乘坐不同等級座位。通過相關變數繪製箱線圖

#選取非空數據繪製箱線圖nembark_fare <- full %>%n filter(PassengerId != 62 & PassengerId != 830)nn# 可視化票價，港口，及乘客客場等級圖nggplot(embark_fare, aes(x = Embarked, y = Fare, fill = factor(Pclass))) +n geom_boxplot() +n geom_hline(aes(yintercept=80), n colour=red, linetype=dashed, lwd=2) +n scale_y_continuous(labels=dollar_format()) +n theme_few()n

由上圖可得，1等因子，票價在80，最接近的是C登船口，可以將缺失值替換為C。

#將C登船口賦值給空值nfull$Embarked[c(62, 830)] <- Cn

2,Fare變數缺失值

#處理Fare變數缺失值na<-which(full$Fare %in% NA)n#輸出缺失值信息nfull[full$PassengerId==a,]n#繪製概率密度圖nggplot(full[full$Pclass == 3 & full$Embarked == S, ], n aes(x = Fare)) +n geom_density(fill = #99d6ff, alpha=0.4) + n geom_vline(aes(xintercept=median(Fare, na.rm=T)),n colour=red, linetype=dashed, lwd=1) +n scale_x_continuous(labels=dollar_format()) +n theme_few()n

此觀測值是由S登船口上船的，第3類客戶，繪製概率密度圖，數據大部分聚集在0-20美金區間，由圖可將中位數替換缺失值。

full$Fare[1044] <- median(full[full$Pclass == 3 &n full$Embarked == S, ]$Fare, na.rm = TRUE)n

3,AGE變數缺失值

# mice包主要用於插補數據，先將變數因子化處理nfactor_vars <- c(PassengerId,Pclass,Sex,Embarked,n Title,Surname,Family,FsizeD)nnfull[factor_vars] <- lapply(full[factor_vars], function(x) as.factor(x))nn# 設置隨機變數種子nset.seed(129)nn# 刪除除創建的變數之外的其他變數，nmethod=隨機森林（rf）nmice_mod <- mice(full[, !names(full) %in% nc(PassengerId,Name,Ticket,Cabin,Family,Surname,Survived)], method=rf) n# 保存完整的輸出nmice_output <- complete(mice_mod)n# 繪製圖形npar(mfrow=c(1,2))nhist(full$Age, freq=F, main=Age: Original Data, n col=darkgreen, ylim=c(0,0.04))nhist(mice_output$Age, freq=F, main=Age: MICE Output, n col=lightgreen, ylim=c(0,0.04))n

# 將模型輸出值替換原數據集年齡的值nfull$Age <- mice_output$Agen

再次觀察下缺失值是否都填補完整。

> sapply(full,function(x) sum(is.na(x)))nPassengerId Survived Pclass Name Sex Age n 0 418 0 0 0 0 n SibSp Parch Ticket Fare Cabin Embarked n 0 0 0 0 1014 0 n Title Surname Fsize Family FsizeD Deck n 0 0 0 0 0 1014 n Child Mother n 0 0 n

現在我們擁有一個完整的數據集。接下來再思考一下，哪些屬性會影響生存狀態。會被優先救助？

ggplot(full[1:891,], aes(Age, fill = factor(Survived))) + n geom_histogram() + n facet_grid(.~Sex) + n theme_few()n

明顯性別變數不同生存率不等，進一步看下具體差異

> table(full$Sex, full$Survived)n n 0 1n female 81 233n male 468 109n> mytable<-xtabs(~Sex+Survived,data=full)nn#頻數列聯表（百分比）n> addmargins(prop.table(mytable)*100)n SurvivednSex 0 1 Sumn female 9.090909 26.150393 35.241302n male 52.525253 12.233446 64.758698n Sum 61.616162 38.383838 100.000000n

船上有38%的人員生存下來，而女性的比例明顯比男性高出一倍。那麼如果是18歲以下的成員呢？

> full$Child[full$Age < 18] <- Childn> full$Child[full$Age >= 18] <- Adultn> table(full$Child, full$Survived)n n 0 1n Adult 481 272n Child 68 70n

nn18歲以下人員獲得了50%的生存機會。再觀察下如果是一位母親，她的生存幾率會不會更高

# 添加變數母親nfull$Mother <- Not Mothernfull$Mother[full$Sex == female & full$Parch > 0 & n full$Age > 18 & full$Title != Miss] <- Mothern> # 母親屬性與生存的列聯表n> table(full$Mother, full$Survived)n n 0 1n Mother 16 39n Not Mother 533 303n

數據表現我們的想法是有理可依的，可添加此兩個變數作為特徵值

#因子化nfull$Child <- factor(full$Child)nfull$Mother <- factor(full$Mother)n

建模：

# 拆分數據集，進行預測訓練ntrain <- full[1:891,]ntest <- full[892:1309,]nn# 設置隨機種子nset.seed(754)nn# 構建模型，保留有用變數nrf_model <- randomForest(factor(Survived) ~ Pclass + Sex + n Age + SibSp + Parch + n Fare + Embarked + Title + n FsizeD+ Child + Mother,n data = train)nn# 使用測試數據集進行預測nprediction <- predict(rf_model, test)nn# 保存測試ID和測試結果生存與否nsolution <- data.frame(PassengerID = test$PassengerId,n Survived = prediction)nn# 寫出結果文件nwrite.csv(solution, file = rf_mod_Solution.csv, row.names = F)n

上傳。可看到排名為3549

在建立模型之前都是在整理數據，整好數據之後可以幫助我們快速的使用更好的方式建立模型。顯然3549的排名並不是非常理想。嘗試使用其他的方式

#利用bagging技術進行建模nBAG=bagging((Survived) ~ Pclass+Age+Sex+SibSp+Parch + n Fare + Embarked + Title + Fsize + FsizeD,n data = train,mfinal=900000)n#對於預測數據做回歸預測nPredictio <- predict(BAG, test)nnsub<- data.frame(PassengerId = test$PassengerId, Survive = Predictio)nsub$Survived[sub$Survive<0.5]<-0nsub$Survived[sub$Survive>=0.5]<-1nsub1=sub[,-2]nsummary(sub1)nn PassengerId Survived n 892 : 1 Min. :0.0000 n 893 : 1 1st Qu.:0.0000 n 894 : 1 Median :0.0000 n 895 : 1 Mean :0.3756 n 896 : 1 3rd Qu.:1.0000 n 897 : 1 Max. :1.0000 n (Other):412 nwrite.csv(sub1, file = "best.csv", row.names = FALSE)n

訓練集有38%的乘客獲得生存機會，此數據已經非常接近了。

上傳結果排名1080。機器學習理論知識，是在此次訓練當中欠缺的地方，之後的學習方向會往這方面偏移。數據分析的思路非常重要，生活經驗有時會誤導我們判斷。比如小孩和母親這兩個特徵值。它們的重要性檢測水平甚至低於登船的港口號。以下是觀測重要性代碼。

#觀察變數重要性nimportance <- importance(rf_model)nvarImportance <- data.frame(Variables = row.names(importance), n Importance = round(n importance[ ,MeanDecreaseGini],2))nn# 對重要性排序nrankImportance <- varImportance %>%n mutate(Rank = paste0(#,dense_rank(desc(Importance))))nn# 畫出重要程度圖nggplot(rankImportance, aes(x = reorder(Variables, Importance), n y = Importance, fill = Importance)) +n geom_bar(stat=identity) + n geom_text(aes(x = Variables, y = 0.5, label = Rank),n hjust=0, vjust=0.55, size = 4, colour = red) +n labs(x = Variables) +n coord_flip() + n theme_few()n

推薦閱讀：