@TOC

1 应用场景

逻辑回归被广泛应用在目标变量是二值变量的场合 (0, 1)

2 公式

– P(y=1|x) 表示 y = 1的概率
– 从而得到 y = 1 对 y = 0 概率的比值

3 模型估计

– 极大似然估计

模型阐释/评估
– 一个解释变量的阐释图

4 读取与整理数据

整理数据格式，将数值变量和因子变量分类

data0=read.csv("accepts.csv",na.strings = "",colClasses = newcolclass)
str(data0)
nrow(data0)
summary(data0)
#定义目标角色
target0="bad_ind"
colclass=sapply(data0,class)
newcolclass=ifelse(colclass=="character","factor","numeric")

#找出data0中变量时numeric的名称
numerics=names(which(sapply(data0,is.numeric)))
length(numerics)

#找出data0是factor的变量名称,集合的差
factors=setdiff(names(which(sapply(data0,is.factor))),target0)
length(factors)

5 变量筛选

5.1集中度和分类数筛选

###取值集中在某一取值
maxper=sapply(data0, function(x) max(table(x)/length(x)))
summary(maxper)

centors=names(which(maxper>0.95))
#分类太多的可以考虑合并
maxperf=sapply(data0[,factors], function(x) max(table(x)/length(x)))
summary(maxperf)

5.2 缺失值分析

nanumber=sapply(data0,function(x) sum(is.na(x)))
summary(nanumber)

#查找大于80%的缺失值
missings3=names(which(nanumber>=0.7*nrow(data0)))

#查找大于50%的缺失值
missings2=names(which(nanumber>=0.5*nrow(data0)))

#查找大于30%的缺失值
missings1=names(which(nanumber>=0.3*nrow(data0)))

修改缺失值的值，数值型用中位数填补，因子型用众数填补

#修改数值型和因子型列缺失值
modifyna=function(x){
  for (i in 1:length(x)) {
    if(sum(is.na(data0[,i]>0))){
      if(class(data0[,i])=="numeric"){
      x[,i][is.na(x[,i])]=median(x[,i],na.rm = T)
      }else{
        a=table(x[,i])
        b=which(a==max(a))
        x[,i][is.na(x[,i])]=names(b)
      }
    }
  }
  x
}
data1=modifyna(data0)

填补后缺失值为0

#将y放在第一列
name=colnames(data1)
match("bad_ind",name)
data2=data1[,c(3,1:2,4:length(data1))]

5.3 频率筛选

#根据业务删除与y无关的

nrow(data2)
colname2=colnames(data2)
match(c("application_id","account_number"),colname2)
data3=data2[,c(-2,-3)]

对过多种类的分类变量进行处理

lvls=sapply(data2[factors],function(x) length(levels(x)))
summary(lvls)


manylvls=names(which(lvls>8))
manylvls
data2$vehicle_make
factors=setdiff(factors,manylvls)

*删除数值类型少的numeric变量，并把二分类变量变为logic因子变量

concret=names(which(lev!=2 & lev<30))
concret
#logic变量
logic=names(which(lev==2))
logic
data3[,logic]=as.data.frame(lapply(data3[,logic],as.factor))
colnames(data3)[2]="vehicle_year"

#确定第一次删除的变量

ignores=c(centors,missings3,manylvls,concret,target0)
input=setdiff(names(data3),ignores)

IV（信息价值）通常应用于评分模型，是基于WOE（证据权重）计算的，
反应的是单个自变量自身与因变量之间的关系，衡量自变量的预测能力，
是筛选变量参考的指标。

计算变量的IV值，ln(bad/good)*(bad-good)

x=data3$age_oldest_tr
target=data3[,1]
summary(x)

x=binning(x,bins=8,method="quantile")
temp1=table(x,target)
temp1

#频数统计表转化为数据框，矩阵的长度为水平的个数
temp2=as.data.frame(matrix(temp1,length(unique(x)),2))
length(unique(temp1))
#计算比例
temp3=as.data.frame(sapply(temp2, function(x) x/sum(x)))

#权重
woe=log(temp3[,2]/temp3[,1])*100
woe

#IV值
iv=sum((temp3[,2]-temp3[,1])*woe)

###自编函数
var.iv=function(x,target){
  if(is.numeric(x)){
    #分箱
    x=binning(x,bins=8,method = "quantile")
    #频数表
    x=table(x,target)
    #转列表
    x=as.data.frame(matrix(x,nrow(x),ncol(x)))
    #比例矩阵
   x=sapply(x,function(x) x/sum(x))
    if(!is.matrix(x)){Iv=0} else{
      #权重
      woe=log(x[,2]/x[,1])*100
      #IV值
      Iv=sum((x[,2]-x[,1])[!is.infinite(woe)]*woe[!is.infinite(woe)])
    }
   return(Iv)
  }
}

var.iv(data3[,12],target)

data4=data3[input]
class(data4)
#计算输入数据的IV值
iv.value=sapply(data4,function(x) var.iv(x,target))
IV=unlist(iv.value)
#增加IV等于0的为因子变量
factors=append(factors,"used_ind")
class(IV)
summary(IV)

#列表转向量
as.vector(unlist(iv.value))

#修改数据框某一列的数据类型
data4$used_ind=as.factor(data4$used_ind)

#选择IV大于10，bad和good差较大的变量
ninput=names(which(IV>10))
length(ninput)

对于分类型变量，也应该筛选IV>10比较合适，这里是为了后面逐步回归比较AIC值，选择大于1的IV值

fvar.iv=function(x,target){
  if(is.factor(x)){
    #分箱
    #频数表
    x=table(x,target)
    #转列表
    x=as.data.frame(matrix(x,nrow(x),ncol(x)))
    #比例矩阵
    x=sapply(x,function(x) x/sum(x))
    
      #权重
      woe=log(x[,2]/x[,1])*100
      #IV值
      fIv=sum((x[,2]-x[,1])[!is.infinite(woe)]*woe[!is.infinite(woe)])
      return(fIv)
    
  }
}
fvar.iv(data4[,1],target)

fiv.value=sapply(data4,function(x) fvar.iv(x,target))
Fiv=unlist(fiv.value)
class(Fiv)

###确定最终进入模型的变量

finput=names(which(Fiv>1))
input=append(ninput,finput)
data_input=data4[,input]
str(data_input)
summary(data_input)

###确定最终进入模型的变量

finput=names(which(Fiv>1))
input=append(ninput,finput)
data_input=data4[,input]
str(data_input)
summary(data_input)

#分割数据设置，训练集和测试集

set.seed(1)
n=nrow(data_input)
train_number=sample(1:n,round(n*0.8))

##构建训练集

train=data_input[train_number,]
train_y=target[train_number]

table(train_y)

test=data_input[-train_number,]
test_y=target[-train_number]

构建模型

logitfit=glm(train_y~.,data=train,family="binomial")
summary(logitfit)

coef(logitfit)

#系数的含义优势比的增加量
sapply(coef(logitfit), function(x) round(exp(x),3))

logitfitstep=step(logitfit)

logitfit=glm(train_y ~ age_oldest_tr + tot_rev_line + fico_score + ltv,data=train,family = "binomial")

通过逐步回归，选择能使AIC最小的变量进入模式

AIC,即赤池信息量准则，是衡量统计模型拟合优良性的一种标准。

●赤池信息量准则建立在熵的概念基础上，可以权衡所估计模型的复杂度和此模型拟合数据的优良性。

AIC= =2k- 2Ln(L)

其中k是模型参数个数，L是似然函数。通常选择AIC最小的模型为最优模型。

●当两个模型之间存在较大差异时，差异主要体现在似然函数项，当似然函数差异不显著时，上式第一项，即模型复杂度则起作用，从而参数个数少的模型是较好的选择。

BIC,即贝叶斯信息准则，与AIC相似，用于模型选择。

训练模型时，增加参数数量,也就是增加模型复杂度,会增大似然函数,但是也会导致过拟合现象,针对该问题，AIC和BIC均引入了与模型参数个数相关的惩罚项，BIC的惩罚项比AIC的大，考虑了样本数量，样本数量过多时，可有效防止模型精度过高造成的模型复杂度过高。

BIC= kLn(n)- 2Ln(L)

… 其中，k为模型参数个数，n为样本数量，L为似然函数。kln(n)惩罚项在维数过大且训练样本数据相对较少的情况下，可以有效避免出现维度灾难现象。

6 对模型进行样本内外预测

library(gtools)
library(gmodels)

#样本内预测

logit_train_fit=fitted(logitfit)
summary(logit_train_fit)
summary(target)

###########################################样本外预测

logit_test_fit=predict(logitfit,newdata = test,type="response") #加responce输出是P，否则是y值
summary(logit_test_fit)

###########################################样本外预测

logit_test_fit=predict(logitfit,newdata = test,type="response") #加responce输出是P，否则是y值
summary(logit_test_fit)

#预测效果检验

##########################################逻辑分类-样本内
p=table(train_y)/length(train_y)

##以p[2]作为临界值划分
logit.cluster=factor(logit_train_fit>p[2],levels=c(F,T),labels=c(0,1))
table(logit.cluster)

#预测准确率
per.pre.train=table(train_y,logit.cluster,dnn=c("Actual","predicted"))
per.pre.train
p11=per.pre.train[2,2]/sum(per.pre.train[2,])
p11

p00=per.pre.train[1,1]/sum(per.pre.train[1,])
p00

accuracy=(per.pre.train[1,1]+per.pre.train[2,2])/sum(per.pre.train)
accuracy

分别查看11，00和总体的正确率

#######################################样本外
#p.pred=exp(logit_test_fit)/(1+exp(logit_test_fit)) 没有responce的概率处理
p.pred=logit_test_fit
summary(p.pred)

pt=table(test_y)/length(test_y)
pt
log(pt[2]/pt[1])
logit.pred=factor(p.pred>p[2],levels=c(F,T),labels=c("0","1"))
table(logit.pred)

#对比预测与真实值

table_pre.test=table(test_y,logit.pred,dnn = c("Actual","predicted"))
table_pre.test

pt11=table_pre.test[2,2]/sum(table_pre.test[2,])
pt11
pt00=table_pre.test[1,1]/sum(table_pre.test[1,])
pt00

accuracyp=(table_pre.test[1,1]+table_pre.test[2,2])/sum(table_pre.test)
accuracyp

这里的预测阈值采用的是训练集中1发生概率的阈值做比较，效果一般

7ROC曲线检验模型

#####################################预测效果检验
library(colorspace)
library(pROC)

modelroc=roc(test_y,logit_test_fit,plot=T)

plot(modelroc,print.auc=T,auc.polygon=T,grid=c(0.1,0.2),grid.col=c("green","red"),max.auc.polygon=T,auc.polygon.col="skyblue",print.thres=T)

modelauc=auc(modelroc)
modelauc

通过roc曲线可以选择最适合的阈值