逻辑回归

本节内容：

1：什么时候用逻辑回归

2：对于脱敏数据，可是自变量又非常多怎么处理-》变量的处理

3：怎么求数据集中的IV值--》IV衡量变量的预测能力

4：创建回归模型

当一系列连续型或类别型预测变量来预测二值型结果变量时，采用Logistic回归是一个非常好用的工具

对逻辑回归参数系数的解释：

在Logistic回归中，响应变量是Y=1的对数优势比（log）。回归系数的含义是当其他预测变量不变时，一单位预测变量

的变化可能引起的响应变量对数优势比的变量。

什么是优势比：简单的栗子

二、对于脱敏数据，可是自变量又非常多怎么处理

#第一步找到哪些是连续型变量  哪些是分类型变量
numerics = names(which(sapply(mtcars,is.numeric)))
factors = setdiff(names(which(sapply(mtcars,is.factor))),target) 
##target是要研究的因子变量

#第二步 找到数据集中每个变量中的最大集中占比是多少
maxper = sapply(mtcars,function(x)max(table(x)/length(x)))
centors = names(which(maxper>0.95))
#0   1 
#12999  23
#如数据里面有这个这样的table数集，0的占比就是99.9%多了 一般我们就要舍弃掉
#这个变量，一个自变量就2个，0占了99%还怎么变的起来

#第三步 找缺失值的数量
nanumber = sapply(data,function(x)sum(is.na(x)))
missages = names(which(nanumber>=0.7*nrow(data)))

##第四步找到水平因子多的，比如说你分析一个1000多行的数据，里面一个自变量他分成了
#500类别，对于这种的我们一般直接舍弃这类变量

lvls = sapply(data[factors],function(x)length(levels(x)))
manylvls = names(which(lvls>15))

##第五步 找到连续型变量中取值较少的数值型  
#如一个自变量是连续型，数据1000行，unique之后这个数据就自由100个独一的数字
lev = sapply(data[numerics],function(x)length(unique(x)))
concret = names(which(lev!=2 & lev<=100))

##第六步 找到只有2个分类的连续变量，其实2个分类一般就是logical值了
logic = names(which(lev==2))


##确定删除的变量
ignores = c(centors,missages,manylvls,concret)

input = setdiff(names(data),ignores) ##祛除不要的自变量的数据

三、怎么求数据集中的IV值--》IV衡量变量的预测能力　

求IV值的函数.var.iv

var.iv = function(x,target){
  if(is.numeric(x)){
    #如何x是数值型的，则对x分段，分8段，按照分位数分
    x<-binning(x,bins=8,method = "quantile")}
  ##做频数统计表
  temp1 = table(x,target)
  ##频数统计表转换为数据框，矩阵的长度为水平的个数
  temp2 = as.data.frame(matrix(temp1,length(unique(x)),2))
  ##计算比例
  temp3 = sapply(temp2,function(x) x/sum(x))
  
  ##如果有不是矩阵的，则直接默认iv=0
  if(!is.matrix(temp3)){iv=0}else{
    woe<-log(temp3[,2]/temp3[,1])*100
    iv<-sum((temp3[,2]-temp3[,1])[!is.infinite(woe)]*woe[!is.infinite(woe)])
    return(iv)
  }
}

library(AER)
data("Affairs")
qq = Affairs

qq$ynaffair[qq$affairs>0]<-1
qq$ynaffair[qq$affairs==0]<-0

##获得IV值 并筛选IV值
#IV值
iv.value = apply(qq,2,var.iv,qq$ynaffair)

##求iv值的统计量如何变量多我们根据IV值进行删除当IV值<12时候我们把变量删除
summary(iv.value)

ivinput = names(which(iv.value>12))
#
#> ivinput
#[1] "age"           "yearsmarried"  "religiousness" "rating" 

四、创建逻辑回归模型

##本次采用的数据集是AER包中的Affaits对于判断是否婚外情的逻辑回归案例
library(AER)
data = Affairs


# 1：了解数据和清洗数据 -------------------------------------------------------------

summary(data)
names(data)
table(data$affairs) ##研究的是自变量婚外情的次数

##对数据婚外次数添加新变量，1表示有婚外情，0表示没婚外情
data$ynaffairs[data$affairs>0] = 1
data$ynaffairs[data$affairs<=0] = 0
data$ynaffairs = factor(data$ynaffairs,levels = c(0,1),labels = c("NO","YES"))

data = data[,-1]  ##去掉原本的affairs变量
names(data)

##这里如果数据是脱敏自变量有非常多是时候，对自变量要进行删除



# 2：分训练集合和测试集 -------------------------------------------------------------
#第一种
train_number= sample(1:nrow(data),round(nrow(data)*0.8))
train = data[train_number,]
test = data[-train_number,]

#第二种
train_number2= sample(2,nrow(data),replace = T,prob = c(0.8,0.1))
train_2 = qq[train_number2,]
test_2 = qq[-train_number2,]


# 3 创建逻辑回归 ----------------------------------------------------------------

fit = glm(train$ynaffairs~.,data=train,family = "binomial")
summary(fit)  ##可以看到有好几个不显著的 去掉这些变量重新拟合

step(fit)  ##可以通过逐步回归

fit = glm(train$ynaffairs ~ gender + age + yearsmarried + religiousness + 
            rating,data=train,family = "binomial")

summary(fit)


# 4：评价这个模型好不好 -------------------------------------------------------------

#通过样本内预测 和 样本外预测
##通过以下的在训练集中没有婚外情占的比例和有婚外情占的比例确定临界点 
##模型的预测值是概率，我们看出有婚外情的是小概率事件
#那么在预测值中p>0.25 标记为0 否则标记为1
p = table(train$ynaffairs)/length(train$ynaffairs)  
#0.7525988 0.2474012 


#样本内预测值并通过概率事件标识0或1
log_fitted = fitted(fit)
fitted_aff = factor(log_fitted>p[2],levels = c(F,T),labels = c(0,1))
table(fitted_aff)

#样本内预测和真实值
apt = table(train$ynaffairs,fitted_aff,dnn=c("actual","predected"))

#         predected
#actual   0   1
# NO     247 115
# YES     44  75

##真实预测1的
apt[2,2]/(apt[2,1]+apt[2,2]) # 0.6302521

##真实预测0的
apt[1,1]/(apt[1,1]+apt[1,2])  #0.6823204

##总的模型预测率
(apt[1,1]+apt[2,2])/sum(apt) ##0.6694387


##通过样本外预测
p.pred = predict(fit,newdata=test,type="response") 
##添加type之后才可以能到预测值的p值
#样本内预测值并通过概率事件标识0或1

pred_p = factor(p.pred>p[2],levels = c(F,T),labels = c("NO","YES"))
table(pred_p)

table(test$ynaffairs)

pre_t= table(pred_p,test$ynaffairs)

#
(pre_t[1,1]+pre_t[2,2])/sum(pre_t)#0.55


##总结对样本内预测达到了67% 对样本外预测达到了55%