用不同方法在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
模型评估的另一种方式：
混淆矩阵相关知识：
预测准确率=正确预测的正反例数 / 总数
误分类率 = 错误预测的正反例数 / 总数
真阳性率(召回率)：正确预测到的正例的数量/实际上正例的个数，在实际应用中覆盖率是很重要的指标，
阳性预测值=正确预测到的正例数 / 预测正例总数
真阴性率 = 正确预测到的负例个数 / 实际负例总数
阴性预测值 = 正确预测到的负例个数 / 预测负例总数
