项目01_Titanic R代码

###############################Step 1 start ##########################

#加载数据源文件
#将空字符，包含空格的字符，NA字符统一处理成缺失值NA
#初始，不将string转换成Factor,影响相关性分析
train <- read.csv("train.csv",na.strings = c(""," ","NA"),stringsAsFactors = FALSE)  
test <-  read.csv("test.csv",na.strings = c(""," ","NA"),stringsAsFactors = FALSE)

#对比train数据集，仅相差一列Survived，对应结果列
str(train)  #891 obs. of  12 variables
str(test)   #418 obs. of  11 variables

#创造test中Survived空列，利用rbind整合成一个大数据集
#一次性加工（缺失值，属性加工，创建新特征等等）训练集数据，测试集
test$Survived <- NA #新增列为逻辑值列
comb <- rbind(train,test) #rbind函数不要求列排序，只需要都存在即可
str(comb) #1309 obs. of  12 variables,合并后Survived跟随训练集的Factor类型
comb$Survived #前891行属于train数据集，后418行属于test数据集

###############################Step 1 stop##########################

###################Step 2 start ###########################
#探索缺失值

#探索数值型变量
#主动忽略Survived列，构造出来的NA值
library("lattice")
library("mice")
md.pattern(comb[,-2]) #输出各列缺失值统计频数表格,仅对NA缺失值统计有效

library("VIM")
aggr(comb[,-2],prop=FALSE,numbers=TRUE) #数值以个数显示
#结果解读，共计缺失1280个值
#Fare列，缺1个值
#Embarked列，缺2个值
#Age列，缺(240+23=263)263个值
#Cabin列，缺(773+240+1=1014)1014个值
aggr(comb[,-2],prop=TRUE,numbers=TRUE) #数值以比例显示
#结果解读，共计缺失1280个值
#Fare列，缺失比例：0.076%
#Embarked列，缺失比例：0.153%
#Age列，缺失比例：(1.757%+18.335%=20.092%)20.092%
#Cabin列，缺失比例：(0.076%+18.335%+59.053%=77.464%)77.464%

##########################step 2 stop ##########################

##########################Step 3 start ########################

#3.1 描述性分析
#对于数值列，五数分析，包含缺失值个数统计
#对于字符列，仅显示观测值个数，列类型为character
summary(comb)
#数值列包含：7列，PassengerId，Survived，Pclass，Age，SibSp，Parch，Fare
#字符列包含：5列，Name，Sex，Ticket，Cabin，Embarked

#3.2 图形探索分析
library(ggplot2)
#3.2.1 Age列对Survived列的相关性
ggplot(train,aes(x=Age,y=PassengerId,color=as.factor(Survived)))+
    geom_point()+
    facet_grid(Sex~.)+
    labs(title="Corr between Age & Survived",x="Age",y="ID",fill="1")+
    scale_color_manual(values=c("blue","red"))+
    theme(legend.title = element_blank(),
          plot.title = element_text(hjust = 0.5)#当R版本下标题无法居中
          )
#图形解读
#0-15岁的间隔中，男，女的存活比率大致相同
#50岁以上，男性死亡比例非常高，女性存活率比例比较高
#对于男性群体而言，50岁以上的死亡率小于20-50年龄段的死亡率

#3.2.2 Pclass列，Embarked列对于Survived列的相关性
ggplot(train[!is.na(train$Embarked),],aes(x=Embarked,y=PassengerId))+
    geom_tile(aes(fill = as.factor(Survived)))+
    facet_grid(.~Pclass)+
    labs(title="Pclass Embarked vs Survived",
         x="Embarked",y="ID")+
    theme( legend.title = element_blank(),
           plot.title = element_text(hjust=0.5)
         )
#图形解读
#从Pclass看，Pclass=3, S港口登录，死亡率最高

ggplot(train[!is.na(train$Embarked),],aes(x=Embarked,y=PassengerId))+
    geom_tile(aes(fill = as.factor(Survived)))+
    facet_grid(.~Sex)+
    labs(title="Sex Embarked vs Survived",
         x="Embarked",y="ID")+
    theme( legend.title = element_blank(),
           plot.title = element_text(hjust=0.5)
    )
#图形解读
#从Sex看，无论男性，还是女性都是S港口登录，死亡率最高

#3.2.3 查看Sex，Age，SibSp，Parch与Survived列之间的相关性
library("grid")
library("vcd")
mosaic(~ Sex + (Age > 15) + (SibSp + Parch > 0) + Survived, 
       data = train[complete.cases(train),],
       shade=TRUE, legend=TRUE)

#图形解读
#马赛克图 mosaic plot
#嵌套矩形面积正比于单元格频率，其中该频率即多维列联表中的频率
#颜色或阴影可表示拟合模型的残差值
#两个以上类别型变量
#皮尔逊残差大于等于2，说明考虑的变量之间存在相互依赖关系。

#sex,Age>15,SibSp+Parch>0,Survived，对应四个维度，相当于多维列表，
#在sex,Age>15,SibSp+Parch>0,三个维度下，分布统计Survived的频数，用不同的残差
#进行图表表现
#考虑Sex,Age>15,SibSp+Parch>0与Survived列之间的相关性

#3.3  探索相关关系
#各自变量(对应其余11列中7个数值型变量)与因变量(Survived)之间的相关关系

#3.3.1 包含缺失值的探索
#corrgram支持类别型，数值型相关性探索,character,Factor列不参与探索
library(corrgram)
corrgram(train,order=TRUE,lower.panel=panel.shade,
         upper.panel=panel.pie,text.panel = panel.txt,
         main="Corrgram of comb intercorrelations"
         
)
#结果解读
#以Survived为直角的，列，行，对应Survived与各列的相关系数
#顺时针，表示正相关，图形扇面越大，相关系数越大
#逆时针，表示负相关，图形扇面越大，相关系数越大
#相关性绝对值排序，扇面大小排序：
#Pclass > Fare > Parch > Age  > Passengrid > SibSp
S_cor <- cor(train[,2],train[,c(1,3,6,7,8,10)],use = "complete.obs")
class(S_cor) #返回matrix
#利用扇形图，来确定相关性大小
library(plotrix)
slices <- abs(S_cor)
lbls   <- colnames(slices)
fan.plot(slices,labels=lbls,main="fan plot") #不加标题，扇形上界超出画布
#确实：Pclass > Fare > Parch > Age  > Passengrid > SibSp
 
#3.3.2 相关系数+P值
#以Survived为直角，上方，向上的直角，图形颜色深浅判断相关性高低
#下方，向下的直角，p的绝对值越大，相关性越高
S_cor_matrix <-  cor(train[,c(1,2,3,6,7,8,10)],use = "complete.obs")
S_cor_matrix_mtest <- cor.mtest(train[,c(1,2,3,6,7,8,10)], conf.level = 0.95,na.action="na.omit")
library(corrplot)
corrplot(S_cor_matrix, method="ellipse",
         p.mat = S_cor_matrix_mtest$p, sig.level = 0.2,
         order = "AOE", 
         type = "upper",
         tl.pos = "d"
)

corrplot(S_cor_matrix, add = TRUE, 
         p.mat = comb_mtest$p, sig.level = 0.2,
         type = "lower", method = "number", 
         order = "AOE",
         diag = FALSE, 
         tl.pos = "n", 
         cl.pos = "n")
#图形解读
#打X表示，P值不显著，近似等于0，可以忽略
#显示相关性排名的末两位，Passengrid ， SibSp，P检验未通过
########################step 3 stop###########################################

########################Step 4 start##########################################
#4.1 特征改造

#4.1.1 Name(包含名，性别/身份.姓),集合Sex,Age，完成分组改造
#针对名字列，完成字符切割，每一个Name分割成3列，第一列为名，第二列为性别/身份，第三列为姓

comb$Title <- sapply(comb$Name,FUN=function(x){strsplit(x,split='[,.]')[[1]][2]})
library("stringr")
comb$Title <- str_trim(comb$Title) #剔除多余的空格
unique(comb$Title) #返回18个值

#基于性别/身份特征，生成对应频数统计表格
T_df    <- as.data.frame(table(comb$Title))
colnames(T_df) <- c("Title","Counts")
T_df <- T_df[order(-T_df$Counts),]
rownames(T_df) <- c(1:18)

#生成辅助表格mid_df
#性别/身份，对应Title字段，人员分组依据
#频数，对应Counts字段，人员分组依据
#Age缺失频数，对应Age_missing字段，
#Age_min,当前组中，年龄最小值
#Age_max,当前组中，年龄最大值

mid_df <- as.data.frame(
          cbind(
              T_df,
              "Age_missing" = sapply(T_df$Title,FUN=function(x) {nrow(comb[comb$Title==x & is.na(comb$Age),])}),
              "Age_min"     = sapply(T_df$Title,FUN=function(x) {min(comb$Age[comb$Title==x],na.rm=TRUE)}),
              "Age_max"     = sapply(T_df$Title,FUN=function(x) {max(comb$Age[comb$Title==x],na.rm=TRUE)})
          )
)
str(mid_df)
mid_df$Title <- as.vector(mid_df$Title) #因子类型处理
mid_df$Title <- str_trim(mid_df$Title)  #空格处理

#   Mr:     For men above 14.5 years
#   Miss:   For girls below and equal to 14.5 years
#   Mrs:    For married women above 14.5 years
#   Master: For boys below and equal to 14.5 years
#   Ms:     For women above 14.5 years, maybe unmarried 综合婚姻基于Miss与Mrs分割出的一组

#精确性分组，修正分裂提取的标记值
comb$Title[comb$Title=="Mr" & !is.na(comb$Age) & comb$Age <=14.5] <- "Master"

q_table <- table(comb$Title,comb$Sex)
c <- rownames(q_table)[q_table[,1]==0 ]#提取非女性是男性Title
c <- c[-c(6,7)] #向量删除
#"Capt"     "Don"      "Jonkheer" "Sir"  
subset(mid_df,Title %in% c,select=c(Title,Age_min,Age_max)) #年龄均在14.5岁以上
comb$Title[comb$Title %in% c] <- "Mr"

c2 <- rownames(q_table)[q_table[,2]==0 ]#提取非男性是女性Title
c2 <- c2[-c(3,6,7)]
#[1] "Dona"         "Lady"         "Mlle"        
#[4] "Mme"          "the Countess"
c_marri <- c2[c(2,5)]
c_unmarri <- c2[-c(2,5)]

#Lady,夫人，the Countess，伯爵夫人，被隐含的婚姻信息
subset(mid_df,Title %in% c2,select=c(Title,Age_min,Age_max)) 
comb$Title[comb$Title %in% c_marri] <- "Mrs"
comb$Title[comb$Title %in% c_unmarri] <- "Ms"

c3 <- rownames(q_table)[q_table[,1]!=0 & q_table[,2]!=0 ]#提取男女都有计数的Title
comb[comb$Title %in% c3,c("Age","Sex")] 
#7个男性，1个缺失值，其余年龄均大于14.5
#1个女性，年龄大于14.5
comb$Title[comb$Title %in% c3 & comb$Sex=="male"] <- "Mr"
comb$Title[comb$Title %in% c3 & comb$Sex=="female"] <- "Ms"

table(comb$Title)
#Master   Miss     Mr    Mrs     Ms 
# 66    260    777    199      7

comb$Title <- as.factor(comb$Title)

#4.1.2 基于SibSp ,Parch ，构造家庭规模FamilySize
#<=3,小家庭
#>3,大家庭
comb$FamilySize <- ifelse(comb$Parch+comb$SibSp+1 <=3,1,0)

#4.1.3 基于Title，Parch判断是否是母亲
comb$Mother <- ifelse(comb$Title=="Mrs" & comb$Parch >0,1,0)

#4.1.4 基于SibSp，Parch，构造single，判断是否单独出行
comb$Single <- ifelse(comb$Parch+comb$SibSp +1 ==1,1,0)

#4.1.5 基于原始Name列，创建FamilyName
comb$FamilyName <- sapply(comb$Name,FUN=function(x){strsplit(x,split='[,.]')[[1]][1]})
Family.Ticket <- comb[comb$Single==0,c("FamilyName","Ticket")]
Family.Ticket <- Family.Ticket[order(Family.Ticket$FamilyName),]
head(Family.Ticket,n=100)
#基于一个家族，共用同一个家族名称，同一个船票号码
#FamilyName+Ticket末尾后3位数字
comb$FamilyName <- paste(comb$FamilyName,str_sub(comb$Ticket,-3,-1),sep="")
head(comb$FamilyName)

#4.1.6 基于每一个人的生还情况，构造家庭整体的生还情况FamilySurvived
Families <- comb[comb$Parch+comb$SibSp >0,] #以2为家庭最小的单位
summary(Families) #Survived列中包含缺失值
Families_Sur <- aggregate(as.numeric(Families$Survived),
                          by=list("FamilyName"=Families$FamilyName),
                          FUN=sum,
                          na.rm=T) #汇总各家庭生还人数
Families_T_Sur <- Families_Sur$FamilyName[Families_Sur$x>0]
comb$FamilySurvived <- 0
comb[apply(comb,1,function(x){
                              ifelse(x["FamilyName"] %in% Families_T_Sur,TRUE,FALSE)
                             }
           ),
     ]$FamilySurvived = 1

#4.1.7 基于Age连续变量,构建类别型变量AgeClass
#NA值暂时忽略，不处理
comb$AgeClass = ifelse(comb$Age<=10,1,
                  ifelse(comb$Age>10 & comb$Age<=20,2,
                         ifelse(comb$Age>20 & comb$Age<=35,3,4))
                  )

##############################step 4 stop###############################

############################step 5 start ################################
#缺失值填充
str(comb)
md.pattern(comb[,-2])
aggr(comb[,-2],prop=FALSE,numbers=TRUE) #数值以个数显示
#读图反馈
#fare列，缺1个
#Embarked列，缺两个
#Age列，缺263个
#Cabin列，缺1014个
#共计缺1280:=1+2+263+1014

#5.1 Embarked列填充
comb[is.na(comb$Embarked),] #同样特征提取，Pclass=1,Survived=1,Cabin=B28,Sex=female
x <- comb$Embarked[!is.na(comb$Survived) &
     comb$Survived==1 & 
     comb$Pclass==1 & 
     comb$Sex=="female" &
     length(grep('^B',comb$Cabin)) > 0
     ]
table(x)
# C  Q  S 
#42  1 46 
#说明S是众数
comb$Embarked[is.na(comb$Embarked)] <- "S"

#5.2 Fare列填充
#决策树填充法
#注意缺失值行其他列不能引入新类型值，否则预测报错
library(rpart)
fit_Fare <- rpart(Fare~Pclass+SibSp+Parch+Age+Embarked+Title,data=comb[!is.na(comb$Fare),],
                  method="anova"
                  )
comb$Fare[is.na(comb$Fare)] <- predict(fit_Fare,newdata=comb[is.na(comb$Fare),])

#5.3 Age列填充
##3.3.2 age列缺失值填充
#缺失263个值，缺失比例为20.092%
#任意值填充容易引起很高的白噪音，降低预测准确率至80%以下

#探索age列与其他列的相关性
library(corrplot)
#计算因变量列，与所有数值型自变量的相关系数
#空值忽略
comb_cor <- cor(comb[,-c(2,4,5,9,11,12,13,14,18)],use = "complete.obs")
#计算相关系数矩阵

#显著性检验，为每一对输入特征产生p值和置信区间
comb_mtest <- cor.mtest(comb[,-c(2,4,5,9,11,12,13,14,18)], conf.level = 0.95,na.action="na.omit")
corrplot(comb_cor, #相关系数矩阵
         method="ellipse",#指定可视化的方法，可以是圆形、方形、椭圆形、数值、阴影、颜色或饼图形
         p.mat = comb_mtest$p,#矩阵的p值，如果空，参数sig.level, insig, pch, pch。上校,pch。cex是无效的。
         sig.level = 0.2,#显著性水平，如果p.mat中的p值大于s .水平，则认为相应的相关系数不显著。
         order = "AOE", #指定相关系数排序的方法，可以是原始顺序(original)、特征向量角序(AOE)、第一主成分顺序(FPC)、层次聚类顺序(hclust)和字母顺序，一般”AOE”排序结果都比”FPC”要好
         type = "upper", #指定展示的方式，可以是完全的、下三角或上三角
         tl.pos = "d" #指定文本标签(变量名称)的位置，当type=full时，默认标签位置在左边和顶部(lt)，当type=lower时，默认标签在左边和对角线(ld)，当type=upper时，默认标签在顶部和对角线，d表示对角线，n表示不添加文本标签
)
#显示结果，检验通过
#SibSp列，相关性最高，颜色最深，abs(p)值最大
#p值显著为0，非常不重要变量，可忽略

#图形探索Age的生存情况
#观察数据集以train数据下
library(ggplot2)
comb_AgeS <- subset(comb[1:891,],!is.na(Age),select=c("Age","Survived"))
ggplot(data=comb_AgeS,aes(x=Age))+
    geom_histogram()+
    facet_wrap(~Survived,nrow=1)

#图形探索pclass分类下，Age有无状况
#Pclass=1,2时，缺失比例比较低
#Pclass=3时，缺失比例比较高
ggplot(data=comb,aes(x=Pclass,fill=!is.na(Age)))+
    geom_bar(position="dodge")+
    labs(title="Pclass vs Age",
         x="Pclass",y="Count",
         fill="Age Null Or Not"
    )

#继续探索age在Pclass=1,2时的生存情况
#观察数据集以train数据下，Pclass=1,2,3
#Pclass=1中大部分乘客都存活下来了
#Pclass=2和Pclass=3中婴儿的存活概率较高
#Pclass=3的乘客最多，但是大部分都没有存活下下来

unique(comb$Pclass) #返回1,2,3
comb_AgeP <- subset(comb[1:891,],Pclass %in% c(1,2,3) & !is.na(Age),select=c("Age","Survived","Pclass"))
summary(comb_AgeP)
ggplot(data=comb_AgeP,aes(x=Age))+
    geom_histogram()+
    facet_wrap(Pclass~Survived,nrow=3)+
    labs(title="Pclass=1,2,Survived, Pclass for age")
#失踪年龄不是随机分布在所有的班级，而是集中在三等乘客中，
#这说明是一个随机失踪(MAR)问题。
#可以使用Mice package使用pmm或预测性均值匹配方法对这些数据进行归算。

comb["Age"]
str(comb)
summary(comb) #目前仅Survived列，Age列,AgeClass存在缺失值
#MI多重插补,针对综合集,排除Survived列
#假如同时存在多列缺失值，mice支持多列同时填充缺失值
imp <- mice(comb[,-c(2,19)],seed=1234)
imp #多重插补，默认形成5个完整的不含缺失值的数据集
nrow(imp$imp$Age) #插补263行
head(imp$imp$Age)
imp_com <- complete(imp,c(seq(1:5))) #多重插补集批量聚合
class(imp_com) #返回data.frame
str(imp_com) #6545 obs. of  17 variables,6545=1309*5
summary(imp_com) #Age列没有任何缺失值

#单独挑选第二次插补集
#Age非空值照抄，空值列pmm方法填充
Age_Filled <- complete(imp,2) 

ggplot(comb,aes(x=Age))+
    geom_density(data=data.frame(comb$PassengerId,Age_Filled),
                                 alpha=0.2,fill="blue")+
    geom_density(data=comb,alpha=0.2,fill="red")+
    labs(title="Age Distribution",x="Age")
#填补后的插补集与原始数据集中的Age分布状况一致
summary(Age_Filled$Age)
comb$Age <- Age_Filled$Age  #利用插补集Age全列，替换comb集Age全列

#5.3.1 根据Age的补充值，同步填充AgeClass列值
comb$AgeClass <-  ifelse(comb$Age<=10,1,
                       ifelse(comb$Age>10 & comb$Age<=20,2,
                              ifelse(comb$Age>20 & comb$Age<=35,3,4))
)

#5.3.2    数据收尾整理

# Title列原值修正，Miss,Ms区间分离
comb$Title[comb$Title=="Miss" & comb$Age >14.5] <- "Ms"
table(comb$Title)
#原始状态：
#Master   Miss     Mr    Mrs     Ms 
# 66    260    777    199      7
#修正后状态：
#Master   Miss     Mr    Mrs     Ms 
#    66     67    777    199    200 
table(comb$Title,comb$Age>14.5)
#        FALSE TRUE
#Master    65    1
#Miss      67    0
#Mr         4  773
#Mrs        3  196
#Ms         0  200

# 生成因子列
comb$Survived <- as.factor(comb$Survived)
comb$Pclass <- as.factor(comb$Pclass)
comb$Sex <- as.factor(comb$Sex)
comb$Embarked <- as.factor(comb$Embarked)
comb$FamilySize <- as.factor(comb$FamilySize)
comb$Mother <- as.factor(comb$Mother)
comb$Single <- as.factor(comb$Single)
comb$FamilyName <- as.factor(comb$FamilyName)
comb$AgeClass <- as.factor(comb$AgeClass)

#5.4  Cabin列值填充，根据来自同一个家庭的其他人的Cabin来进行填充
comb$CabinNo <- sapply(comb$Cabin,function(x){substr(x,1,1)})
comb$CabinNo <- as.factor(comb$CabinNo)
table(is.na(comb$Cabin))
#FALSE  TRUE 
#  295  1014 
#缺失值为1014个

#提取来自同一个FamilyName，对应的CabinNo
Fname_CabinNo <- unique(comb[!is.na(comb$CabinNo)& comb$Parch+comb$SibSp >=1,c("FamilyName","CabinNo")])
str(Fname_CabinNo) #返回83 obs. of  2 variable
head(Fname_CabinNo)

#将跟家庭一同出行，并且缺失CabinNo信息，利用家庭CabinNo进行填充
#显示4个缺失值，可以被填充
Passid <- comb$PassengerId[is.na(comb$CabinNo) & comb$FamilyName %in% Fname_CabinNo$FamilyName]
for (i in 1:length(Passid)){
    myFn <- comb$FamilyName[comb$PassengerId == Passid[i]] 
    myCabin <- Fname_CabinNo[Fname_CabinNo$FamilyName == myFn,"CabinNo"]
    comb$CabinNo[comb$PassengerId == Passid[i]] <- myCabin
}
comb[comb$PassengerId %in% Passid,"CabinNo"]

#剩余1010个Cabin缺失值通过来Pclass值来填充

n_Pclass <- addmargins(table(comb[is.na(comb$CabinNo),c("Pclass")]))
class(n_Pclass) #返回"table" "array"
# 1    2    3   Sum 
#65  254  691  1010 
#缺失值中
#pclass=1,出现65个缺失值
#pclass=2,出现254个缺失值
#pclass=3,出现691个缺失值

nn_Pclass <- addmargins(table(comb$CabinNo,comb$Pclass))
class(nn_Pclass) #返回"table" "array"

#      1   2   3 Sum
#A    22   0   0  22
#B    65   0   0  65
#C    96   0   0  96
#D    40   6   0  46
#E    34   4   3  41
#F     0  13  10  23
#G     0   0   5   5
#T     1   0   0   1
#Sum 258  23  18 299

#结合n_Pclass，nn_Pclass汇总各类出现个数
#Type Pclass=1    Pclass=2    Pclass=3    sum
#A       22    0    0      22
#B       65    0    0      65
#C       96    0    2      98
#D       40    6    0      46
#E       34    4    3      41
#F       0    13    8      21
#G       0    0    5      5
#T       1    0    0      1
#NA       65    254    691     1010
#sum  323    277    709     1309

#计算各类别在Pclass=1,2,3时的出现概率
#            Pclass=1          Pclass=2          Pclass=3
#非空占比    Type/(sum-NA)    Type/(sum-NA)    Type/(sum-NA)
#    A             0.09           0                0
#    B             0.25           0                0
#    C             0.37           0              0.11
#    D             0.16          0                0
#    E             0.13           0              0.17
#    F              0               1              0.44
#    G              0               0              0.28

#利用非空时概率，估计空时的各类大概对应的个数
#空占比    Pclass=1    Pclass=2    Pclass=3
#A           6             0          0
#B          16           0           0
#C        24           0          76
#D          10           0          0
#E           8           0          117
#F           0           254          304
#G           0           0          193
#T           1           0          1
#sum      65           254          691

N_P1 <- c(6,16,24,10,8,0,0,1)
N_P2 <- c(0,0,0,0,0,254,0,0)
N_P3 <- c(0,0,76,0,117,304,193,1)

set.seed(1234)
for(i in 1:nlevels(comb$Pclass)){
    pid <- comb[is.na(comb$CabinNo) & comb$Pclass==i,"PassengerId"]
    #对PassengerId完成随机重排序
    #代替sample函数，多次取，未必能覆盖提取全部数据
    pid <- sample(x=pid) 
    N_P <-switch(i,N_P1,N_P2,N_P3)
    index_start <- 1 #定义pid向量的起始下标
        for (j in 1:nlevels(comb$CabinNo)){
        if(N_P[j]==0){
            next
        }else{
            index_end   <- index_start + N_P[j] -1 #定义pid向量的结束下标
            Pid_pick <- pid[index_start:index_end]
            Cabin_v <- rep(levels(comb$CabinNo)[j],length(Pid_pick))
            comb$CabinNo[comb$PassengerId %in% Pid_pick] <- Cabin_v
            index_start <- index_end + 1
        }
    }
}
nn_Pclass <- addmargins(table(comb$CabinNo,comb$Pclass))
class(nn_Pclass) #返回"table" "array"
#       1    2    3  Sum
#A     28    0    0   28
#B     81    0    0   81
#C    120    0   76  196
#D     50    6    0   56
#E     42    4  120  166
#F      0  267  314  581
#G      0    0  198  198
#T      2    0    1    3
#Sum  323  277  709 1309

comb_bakeup <- comb
#综合数据，保留特征改造，删除特征原始
#Name索引位置4,Ticket索引位置9,Cabin索引位置11
#也可以使用comb$Name <- NULL来删除指定列
comb <- comb[,-c(4,9,11)]
summary(comb)
str(comb)

#拆分训练集，测试集
train <- comb[1:891,]
str(train)
test  <- comb[892:1309,]
str(test)
############################step 5 stop  ################################

############################step 6 start ###############################

####6.1  训练集数据预处理########
train$PassengerId <- NULL #等价于删除列
str(train)

#train_sub_male 做为男性数据集
train_sub_male <- subset(train,Sex=="male")
str(train_sub_male)
#577 obs. of  16 variables
train_sub_male$Sex <- NULL
train_sub_male$Mother <- NULL
train_sub_male$Title <- droplevels(train_sub_male$Title)
#droplevels函数，删除数据集中未使用的levels
#nlevels(trian_sub_male$Title)=2，返回水平数 2
#levels(trian_sub_male$Title) 返回水平名称："Master" "Mr"
#class(trian_sub_male$Title) #返回factor
#有助于提高准确性

#train_sub_female 做为女性数据集
train_sub_female <- subset(train,Sex=="female")
str(train_sub_female)
#314 obs. of  16 variables
train_sub_female$Sex <- NULL #有助于提高准确性
#nlevels(trian_sub_female$Title)=5，返回水平数 5
#levels(trian_sub_female$Title) 返回水平名称："Master" "Miss"   "Mr"     "Mrs"    "Ms"
#class(trian_sub_female$Title) #返回factor
train_sub_female$Title <- droplevels(train_sub_female$Title)
#nlevels(trian_sub_female$Title)=3，返回水平数 3
#levels(trian_sub_female$Title) 返回水平名称："Miss" "Mrs"  "Ms"
#class(trian_sub_female$Title) #返回factor

#性别内部继续衍射测试集，验证集
set.seed(1234)
library("caTools")
split_v <- sample.split(train_sub_male,SplitRatio = 0.75)
class(split_v) #返回，logical
length(split_v) #返回，14

male_train <-  subset(train_sub_male,split_v==TRUE)
str(male_train) #413 obs. of  14 variables
male_test  <-  subset(train_sub_male,split_v==FALSE)
str(male_test)  #164 obs. of  14 variables
#413+164=577,对应trian_sub_male的观察值总量

set.seed(1234)
split_v2 <- sample.split(train_sub_female,SplitRatio = 0.75)
class(split_v2) #返回，logical
length(split_v2) #返回，15
female_train <- train_sub_female[split_v2,]
str(female_train) #230 obs. of  15 variables
female_test <- train_sub_female[!split_v2,]
str(female_test) #84 obs. of  15 variables
#230+84=314，对应train_sub_female的观察值总量

####6.2  训练集上建立模型########

###6.2.1 广义线性模型#######

############################简单logistic 回归#####################
fit_glm <- glm(Survived~.,data=train,family=binomial(link='logit'),control=list(maxit=100))
#1: glm.fit: fitted probabilities numerically 0 or 1 occurred
#所以当样本数据完全可分时，logistic回归往往会导致过拟合的问题，即出现第二个警告：
#拟合概率算出来的概率为0或1。
#出现了第二个警告后的logistic模型进行预测时往往是不适用的，
#对于这种线性可分的样本数据，其实直接使用规则判断的方法则简单且适用（如当pl<2.5时则直接判断为setosa类，pl>2.5时判断为versicolor类）。
summary(fit_glm)
#Number of Fisher Scoring iterations: 21
#仅返回部分结果，并未返回control=list(maxit=100)所选择的迭代次数
#并未返回R方，F检验，P值等选项
AIC(fit_glm) #返回1728.262

library("MASS")
fit_glm_stepAIC <- stepAIC(fit_glm,direction = "both")
AIC(fit_glm_stepAIC) #逐步回归以后，AIC=579.198
summary(fit_glm_stepAIC) #基本与summary(fit_glm)类似
#Call:
#glm(formula = Survived ~ Pclass + Age + Title + FamilySize + 
#        Single + FamilySurvived, family = binomial(link = "logit"), 
#    data = train, control = list(maxit = 100))
#根据AIC 赤池信息准则，
#模型最后选择Pclass，Age，Title，FamilySize，Single，FamilySurvived
#6个自变量参与建模。

#比较模型
#将模型fit_glm_setpAIC嵌套在模型fit_glm中
#anova函数还同时对是否应该添加到线性模型中进行了检验
#根据最优模型选择条件，精度相同时，选模型简洁的，所有不需要将两个变量添加到线性模型中
#P值，F值均未出现
anova(fit_glm,fit_glm_stepAIC)

#检查VIF 方差膨胀因子
#检查多重共线性
library("carData")
library("car")
vif(fit_glm_stepAIC)

pred_fit_glm_stepAIC <- predict(fit_glm_stepAIC,train)
pred_fit_glm_stepAIC[pred_fit_glm_stepAIC <= 0.5] <- 0
pred_fit_glm_setpAIC[pred_fit_glm_stepAIC > 0.5] <- 1
CF_glm_stepAIC <- confusionMatrix(as.factor(train$Survived),as.factor(pred_fit_glm_stepAIC))
logistic_Acc <- CF_glm_stepAIC$overall[1] #正确率 0.8731762 

result <- predict(fit_glm_stepAIC,test)
result[result <= 0.5] <- 0
result[result > 0.5] <- 1
table(result)
############################简单logistic 回归#####################

###6.2.1 广义线性模型#######

####6.2  训练集上建立模型########
#目前最好用的拟合广义线性模型的 R package 是 glmnet
#以下五类模型的变量选择可采用R语言的glmnet包来解决。这五类模型分别是：
#family 规定了回归模型的类型：
#family="gaussian" 适用于一维连续因变量（univariate）
#1. 二分类logistic回归模型

#family="mgaussian" 适用于多维连续因变量（multivariate）
#2. 多分类logistic回归模型

#family="poisson" 适用于非负次数因变量（count）
#3. Possion模型

#family="binomial" 适用于二元离散因变量（binary）
#4. Cox比例风险模型

#family="multinomial" 适用于多元离散因变量（category）
#5. SVM

#这两个都是正则化的手段。
#LASSO是基于回归系数的一范数，Ridge是基于回归系数的二范数的平方。
#模型中有很多变量对模型都有些许影响，那么用Ridge；
#如果你的模型中只有少量变量对模型很大影响，那么用LASSO。
#LASSO可以使得很多变量的系数为0（相当于降维），但是Ridge却不能。
#因为Ridge计算起来更快，所以当数据量特别大的时候，更倾向于用Ridge。
#最万能的方法是用LASSO和Ridge都试一试，比较两者Cross Validation的结果。
#你可以用elastic net, 多了一个超参，但是结合了both ridge and LASSO
#LASSO可以降维，特征选择，可以简化模型；而Ridge不会进行特征选择

#男性生存率预测模型
independent_var_m <- male_train[,2:14] #自变量集
str(independent_var_m)
class(independent_var_m) #data.frame
depentdent_m <- male_train[,1]#因变量集
class(depentdent_m) #返回facotr
#自变量，因变量过拟合问题处理
install.packages("glmnet")
library("glmnet")

#岭回归与Lasso回归的出现是为了解决线性回归出现的过拟合以及
#在通过正规方程方法求解θ的过程中出现的x转置乘以x不可逆这两类问题的
#其中λ称为正则化参数，如果λ选取过大，会把所有参数θ均最小化，造成欠拟合，
#如果λ选取过小，会导致对过拟合问题解决不当，因此λ的选取是一个技术活。 
#岭回归与Lasso回归最大的区别在于岭回归引入的是L2范数惩罚项，
#Lasso回归引入的是L1范数惩罚项，Lasso回归能够使得损失函数中的许多θ均变成0，这点要优于岭回归，
#因为岭回归是要所有的θ均存在的，这样计算量Lasso回归将远远小于岭回归。
#Lasso回归最终会趋于一条直线，原因就在于好多θ值已经均为0，而岭回归却有一定平滑度，因为所有的θ值均存在。

set.seed(1234)
independent_var_m  <- data.matrix(independent_var_m)
glm_fit_male.ridge <- cv.glmnet(independent_var_m,depentdent_m,
                                family="binomial", #对应Logistic 回归
                                alpha=0, #等于0，对应岭回归
                                type.measure="class")
#抛异常
#Error in lognet(x, is.sparse, ix, jx, y, weights, offset, alpha, nobs,  : 
#(list) object cannot be coerced to type 'double
#解决办法：x,必须是matrix类型,y,必须是因子类型
#glmnet 只能接受数值矩阵作为模型输入，如果自变量中有离散变量的话，
#需要把这一列离散变量转化为几列只含有 0 和 1 的向量， 这个过程叫做 One Hot Encoding
print(glm_fit_male.ridge)
#Df       %Dev        Lambda
#[1,] 13 -4.121e-15 160.60000
#[2,] 13  1.047e-03 146.30000
#每一行代表了一个模型。
#列 Df 是自由度，代表了非零的线性模型拟合系数的个数。 
#列 %Dev 代表了由模型解释的残差的比例，对于线性模型来说就是模型拟合的 R^2（R-squred）。 
#它在 0 和 1 之间，越接近 1 说明模型的表现越好，如果是 0，说明模型的预测结果还不如直接把因变量的均值作为预测值来的有效。
#列 Lambda 当然就是每个模型对应的 λ 值。 
#随着 λ 的变小，越来越多的自变量被模型接纳进来，%Dev 也越来越大

glm_fit_male.Lasso <- cv.glmnet(independent_var_m,depentdent_m,
                                family="binomial",
                                alpha=1, #等于1，对应Lasso回归
                                type.measure="class")

#type.measure 是用来指定交叉验证选取模型时希望最小化的目标参量，对于 Logistic 回归有以下几种选择：
#type.measure="deviance" 使用 deviance，即 - 2 倍的 log-likelihood
#type.measure="mse" 使用拟合因变量与实际应变量的 mean squred error
#type.measure="mae" 使用 mean absolute error
#type.measure="class" 使用模型分类的错误率（missclassification error）
#type.measure="auc" 使用 area under ROC curve，是现在最流行的综合考量模型性能的一种参数

p_bakeup <- par(no.readonly = TRUE)
par(mfrow=c(1,2))

plot(glm_fit_male.ridge,main="Ridge")
plot(glm_fit_male.Lasso,main="Lasso")

#lambda.min是指交叉验证中使得均方误差最小的那个值λ
coef(glm_fit_male.ridge,s="lambda.min")
coef(glm_fit_male.Lasso,s="lambda.min")

#Ridge模型给出的错误率更低,最终选用Ridge模型
#在male_train,训练集，查看正确率
pre_glm_fit_male <- predict(glm_fit_male.ridge,independent_var_m,type="class")
table(male_train$Survived,pre_glm_fit_male,male_train$Title)
#male_test，验证集，  查看准确率
pre_glm_fit_male <- predict(glm_fit_male.ridge,data.matrix(male_test[,2:14]),type="class")
table(male_test$Survived,pre_glm_fit_male,male_test$Title)
#男性Ridge模型，在Master类别<=14.5的正确率更高，Mr类别>15的正确率更低

#女性生存率预测模型
str(female_train)
#230 obs. of  15 variables
independent_var_fe <- female_train[,2:15] #自变量集
str(independent_var_fe)
class(independent_var_fe) #data.frame
depentdent_fe <- female_train[,1]#因变量集
class(depentdent_fe) #返回facotr
#自变量，因变量过拟合问题处理
set.seed(1234)
independent_var_fe  <- data.matrix(independent_var_fe)
glm_fit_female.ridge <- cv.glmnet(independent_var_fe,depentdent_fe,
                                  family="binomial", #对应Logistic 回归
                                  alpha=0, #等于0，对应岭回归
                                  type.measure="class")
print(glm_fit_female.ridge)

glm_fit_female.Lasso <- cv.glmnet(independent_var_fe,depentdent_fe,
                                  family="binomial",
                                  alpha=1, #等于1，对应Lasso回归
                                  type.measure="class")

plot(glm_fit_female.ridge,main="Ridge")
plot(glm_fit_female.Lasso,main="Lasso")

#lambda.min是指交叉验证中使得均方误差最小的那个值λ
coef(glm_fit_female.ridge,s="lambda.min")
coef(glm_fit_female.Lasso,s="lambda.min")

#Ridge模型给出的错误率更低,最终选用Ridge模型
#在female_train,训练集，查看正确率
pre_glm_fit_female <- predict(glm_fit_female.ridge,independent_var_fe,type="class")
levels(pre_glm_fit_female) #返回NULL
table(female_train$Survived,pre_glm_fit_female,female_train$Title)

#混淆矩阵展示测试集上正确率
library("caret")
#confusionMatrix(female_train$Survived,pre_glm_fit_female)
#Error: `data` and `reference` should be factors with the same levels.
#强制转成因子型，并保持同样的因子水平
confusionMatrix(female_train$Survived,as.factor(pre_glm_fit_female))

#female_test，验证集，  查看准确率
pre_glm_fit_female <- predict(glm_fit_female.ridge,data.matrix(female_test[,2:15]),type="class")
table(female_test$Survived,pre_glm_fit_female,female_test$Title)
confusionMatrix(female_test$Survived,as.factor(pre_glm_fit_female))
#女性Ridge模型，在Miss类别,Mrs的正确率更高，Ms类别的正确率更低

############################step 6 top   ###############################

############################step 7 start ###############################

####7.1  测试集数据预处理########

#train_sub_male 做为男性数据集
test_sub_male <- subset(test,Sex=="male")
str(test_sub_male)
#266 obs. of  17 variables
test_sub_male$Sex <- NULL
test_sub_male$Mother <- NULL
test_sub_male$Title <- droplevels(test_sub_male$Title)

#train_sub_female 做为女性数据集
test_sub_female <- subset(test,Sex=="female")
str(test_sub_female)
#152 obs. of  17 variables
test_sub_female$Sex <- NULL
test_sub_female$Title <- droplevels(test_sub_female$Title)

#基于男性模型
#在test_sub_male，测试集，预测
#266 obs. of  15 variables
str(test_sub_male)
pre_glm_fit_male <- predict(glm_fit_male.ridge,data.matrix(test_sub_male[,3:15]),type="class")
table(pre_glm_fit_male,test_sub_male$Title)

#基于女性模型
#在test_sub_female，测试集，预测
pre_glm_fit_female <- predict(glm_fit_female.ridge,data.matrix(test_sub_female[,3:16]),type="class")
table(pre_glm_fit_female,test_sub_female$Title)

#统一男性预测值，女性预测值为一个结果集
#注意as.data.frame,以及data.frame的特定使用场合
male_result <- data.frame(PassengerId=test_sub_male$PassengerId,Survived=pre_glm_fit_male)
str(male_result)
female_result <- data.frame(PassengerId=test_sub_female$PassengerId,Survived=pre_glm_fit_female)
str(female_result)
result <- rbind(male_result,female_result) 
str(result)
names(result) <- c("PassengerId","Survived")
result <- result[order(result$PassengerId),]
table(result$Survived)

############################step 7 stop ###############################

############################step 8 start ##############################
#结果集写文件
filepath <- paste("D:\\Script\\R\\kaggle\\Titanic\\gender_submission.csv",sep="")
write.csv(result,filepath,row.names=FALSE,quote=FALSE)

#kaggle目前排名 606/10129,top 6%，最优得分0.81339，
############################step 8 stop ###############################