R语言回归分析实例

本文包含多元线性回归及逻辑回归两种算法，个人实践操作，希望能与大家一起交流分享，如有描述不当之处，欢迎并多谢指正。

#一、公共部分，加载并审核数据、设置数据分区
#-----------------------------
#1、设置工作目录
setwd("E:/R语言建模实例")
#-----------------------------
#2、加载数据，header=TRUE表示读取表头，sep为分隔符设置
userdata <- read.table(file="用户清单.txt",header = TRUE,sep = ",")
#-----------------------------
#3、查看统计量并进行数据审核
vsummary<-summary(userdata)
View(vsummary)

#数据审核相关函数说明
#summary()函数可以获取描述性统计量，可以提供最小值、最大值、四分位数和数值型变量的均值，以及因子向量和逻辑型向量的频数统计
#misc包中的describe()函数，可返回变量和观测的数量、缺失值和唯一值的数目、平均值、分位数，以及五个最大的值和五个最小的值
#psych包中的describe()函数，psych包也拥有一个名为describe()的函数，它可以计算非缺失值的数量、平均数、标准差、中位数、截尾均值、绝对中位差、最小值、最大值、值域、偏度、峰度和平均值的标准误
#pastecs包中的stat.desc()的函数，可以计算种类繁多的描述性统计量。使用格式为：stat.desc（x,basic=TRUE,desc=TRUE,norm=FALSE,p=0.95）
#其中的x是一个数据框或时间序列。若basic=TRUE（默认值），则计算其中所有值、空值、缺失值的数量，以及最小值、最大值、值域，还有总和。若desc=TRUE（同样也是默认值），则计算中位数、平均数、平均数的标准误、平均数置信度为95%的置信区间、方差、标准差以及变异系数。最后，若norm=TRUE（不是默认的），则返回正态分布统计量，包括偏度和峰度（以及它们的统计显著程度）和Shapiro–Wilk正态检验结果
#str()函数，以简洁的方式显示对象的数据结构及内容，可以查看数据框中每个变量的属性
#attributes()函数，可以提取对象除长度和模式以外的各种属性
library(psych)
vdescr<-describe(userdata)
View(vdescr)

#用View查看结果，依次是变量名称、顺序号、元素数量、平均值、标准差、中位数、截尾均值、绝对中位差、最小值、最大值、值域、偏度、峰度和平均值的标准误差
#-----------------------------
#4、相关性矩阵,变量较多，最好是将矩阵保留到变量中，用View()函数查看
vcor<-cor(userdata)
View(vcor)

#用离网标识来举例，可以从结果上查看相关系数比较高的一些指标，
#近三月通话次数波动率为负相关，近三个月欠费次数、当月0通话天数为正相关
#-----------------------------
#5、去掉用户ID，筛选完后可以用str()函数查看一下筛选情况
newusdata<-userdata[,c(-1)]#-1中，-号为排除变量，1表示排除第一列
#6、分箱，分出训练集和测试集
set.seed(1234)
#样本抽取，分两个箱，replace=TRUE为允许放回抽取,FALSE为一次性抽取,prob为抽取概率
ind <- sample(2, nrow(newusdata), replace=TRUE, prob=c(0.5, 0.5))

#获取训练集
trainData <- na.omit(newusdata[ind==1,]) #顺便使用na.omit()函数清洗无效数据
#获取测试集
testData <- na.omit(newusdata[ind==2,]) #顺便使用na.omit()函数清洗无效数据
#可以统计一下训练集数据
nrow(trainData)
#[1] 37528
nrow(testData)
#[1] 37340
#二、多元线性回归
#我们先来看一下ARPU和业务使用量等指标之间的线性回归

#7、训练多元线性回归模型，预测出账收入，先构建ARPU与所有其它变量的线性关系模型
fit_all<-lm(出账收入~., data = trainData)
#8、观察并解读模型统计量
summary(fit_all)
#--------模型fit_all的统计量 start------------
# Call:
#   lm(formula = 出账收入 ~ ., data = trainData)
#
# Residuals:#残差，分别是最小值、第一分位数、中位数、第三分位数、最大值，残差值越小，说明模型拟合度越好，
# 以下残差值的两端极值较大，但从中位数、第一分位数、第三分位数看还是可以
# 理想的情况下，回归残差应该有一个完美的正态分布。1Q和3Q应该有大约相同的幅度，
# 中位数符号则表示数据的倾斜方向，本例中为略向右倾斜
#   Min      1Q  Median      3Q     Max
# -333.54   -2.82    0.81    4.32  182.07
#
# Coefficients:  #系数，分别是系数估计值、标准误差、t值、P-value值（用于T检验显著性标记）
# Estimate是由普通最小二乘法计算出的估计回归系数
# Std.Err为估计系数的标准误差
# t value：t统计量，t检验是用t分布理论来推论差异发生的概率，评价差异是否显著
# Pr是一个概率，它估计系数不显著的可能性，越小越好
#                             Estimate        Std. Error t value Pr(>|t|)  
# (Intercept#截距)                 9.832e-01  3.244e-01    3.031 0.002440 **
# 是否集团客户                  -4.893e-01  1.799e-01   -2.721 0.006520 **
# 是否融合业务                  -1.819e+00  1.545e-01  -11.773  < 2e-16 ***
# 用户当前积分余额              -3.531e-04  3.373e-05  -10.467  < 2e-16 ***
# 终端是否识别                   1.161e-01  1.828e-01    0.636 0.525089   
# 终端类型                       4.762e-01  2.056e+00    0.232 0.816832   
# 是否4G                        -8.803e-02  1.474e-01   -0.597 0.550394   
# 是否智能机                    -9.921e-01  2.051e+00   -0.484 0.628537   
# 入网渠道类型                   5.062e-05  2.850e-05    1.776 0.075787 . 
#......
# 识别出的短号码个数            -8.399e-03  3.327e-02   -0.252 0.800677   
# 合约剩余月数                   1.796e-03  3.835e-05   46.823  < 2e-16 ***
#   是否离网                       5.442e-02  1.807e-01    0.301 0.763357   
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# # Signif. codes 变量的显著性标记，***为最重要，无*则表示变量不可用或不相关
# Residual standard error: 10.23 on 37452 degrees of freedom #表示残差的标准差的自由度
# Multiple R-squared:  0.9745,#为相关系数R^2的检验，表示回归模型所能解释的响应变量的方差比例，越接近1则越显著，越大越好，这里表示解释出账收入的方差为0.9745,剩余0.0255是不能解释的
# Adjusted R-squared:  0.9744 #为相关系数的修正系数，一般多元回归模型建议用此项，它考虑了模型中变量的数目，能实际地评估模型的有效性
# F-statistic: 1.907e+04 on 75 and 37452 DF,  p-value: < 2.2e-16 #表示F统计量，评估模型是否显著，p值小于0.05则表明该模型是显著的，可靠大于0.05则模型不显著，其它的变量也就不用看了
#--------模型fit_all的统计量 end------------
#9、绘制模型散点图
par(mfrow=c(2,2))#2*2的图形分布
plot(fit_all)

#一共四个图：
#1)残差图与拟合图"Residuals vs Fitted，该曲线尽量拟合所有点，若因变量与自变量线性相关，那么残差值与预测（拟合）值就没有任何系统关联
#2)正态Q-Q图"Normal Q-Q"：是在正态分布对应的值下，标准化残差的概率图。若满足正态假设，那么图上的点应该落在呈45度角的直线上；若不是如此，那么就违反了正态性的假设
#3)位置尺度图"Scale-Location"：若满足不变方差假设，那么在位置尺度图中，水平线周围的点应该随机分布。
#4)残差与杠杆图"Residuals vs Leverage"：提供了你可能关注的单个观测点的信息。从图形可以鉴别出离群点、高杠杆值点和强影响点,
#  数字表示异常观测点
# 从所有变量的模型诊断图上来看，效果并不理想，理想的：残差拟合图上的点保持随机分布，正态QQ图上的点应该在45度斜线上，表示满足正态分布
#-----------------------------
# 看一下模型的估计系数置信区间
coef(fit_all)
confint(fit_all)
#                                   2.5 %        97.5 %
#   (Intercept)                    3.473883e-01  1.619087e+00
# 是否集团客户                  -8.418042e-01 -1.367810e-01
# 是否融合业务                  -2.121297e+00 -1.515796e+00
# 用户当前积分余额              -4.191843e-04 -2.869512e-04
# 终端是否识别                  -2.420636e-01  4.743591e-01
# 终端类型                      -3.553723e+00  4.506193e+00
# 是否4G                        -3.769527e-01  2.008956e-01
# ......
# 呼转次数                                 NA            NA#这是有异常值
# 呼转时长                                 NA            NA
# ......
# 出帐收入下降比例.比上月.       7.005021e-04  2.022073e-03
# 本网通话时长占比               6.052787e-01  1.622902e+00
# 连续3月平均上网时长.分.       -5.951849e-05 -3.290378e-05
# 计费原始出账收入               8.050673e-01  8.141825e-01
# ......
# 是否离网                      -2.998516e-01  4.086902e-01

#问题比较严重，在置信区间内，符合相反表示该变量系数可能为0，就是说这变量可以被干掉了，另外NA异常的也可以干掉了
#-----------------------------

# 结合统计量来看，变量过多，需要进行筛选
#10、进行模型变量筛选
#人工筛选？NO，近80个指标筛选太麻烦
full.model<-lm(出账收入~., data = trainData)
reduced.model<-step(full.model,direction = "backward")
# 后向逐步回归法，开始有很多的变量，一步一步收敛，注意观察AIC值的变化
# 收敛过后的指标，进行统计量观察
summary(reduced.model)
# 针对统计量进行观察，进一步收敛模型变量
# tsData<-trainData[c("是否融合业务","用户当前积分余额","是否智能机","合约客户","合约类型","账户余额",
#                 "月末欠费金额","近三个月欠费次数","近三个月缴费金额","赠送费用",
#                 "赠送费用占出账收入比例","总计费时长","被叫计费时长","漫游通话次数",
#                 "主叫通话次数占比","上网计费流量","出帐收入下降比例.比上月.","当月0通话天数",
#                 "连续3月平均出账收入","连续3月平均总通话次数","连续3月平均主叫次数",
#                 "连续3月平均主叫计费时长","连续3月平均被叫次数","连续3月平均点对点短信发送条数",
#                 "连续3月平均上网计费流量","近三月出账收入波动率","是否集团拍照",
#                 "X12个月累计出账收入","当月是否停机","本网通话时长占比",
#                 "连续3月平均上网时长.分.","计费原始出账收入","合约剩余月数")]
# summary(tsData)
bwdmodel<-lm(出账收入~是否融合业务 + 用户当前积分余额 +
                   是否智能机 +  合约客户 +
                   合约类型 + 账户余额 + 月末欠费金额 + 近三个月欠费次数 + 近三个月缴费金额 +
                   赠送费用 + 赠送费用占出账收入比例 + 总计费时长 +
                   被叫计费时长 + 漫游通话次数 + 主叫通话次数占比 + 上网计费流量 +
                   出帐收入下降比例.比上月. + 当月0通话天数 + 连续3月平均出账收入 +
                   连续3月平均总通话次数 + 连续3月平均主叫次数 + 连续3月平均主叫计费时长 +
                   连续3月平均被叫次数 + 连续3月平均点对点短信发送条数 + 连续3月平均上网计费流量 +
                   近三月出账收入波动率 +
                   是否集团拍照 + X12个月累计出账收入 + 当月是否停机 + 本网通话时长占比 +
                   连续3月平均上网时长.分. + 计费原始出账收入 + 合约剩余月数, data = trainData)
summary(bwdmodel)
#这时候发现指标的标记还是存在无*的情况，可酌情进一步筛选
#-----------------------------
# 再看一下收敛后模型的估计系数及其置信区间
coef(bwdmodel)
confint(bwdmodel)
#                                   2.5 %        97.5 %
#   (Intercept)                    4.505271e-01  1.384207e+00
# 是否融合业务                  -2.094359e+00 -1.522474e+00
# 用户当前积分余额              -4.258241e-04 -2.945201e-04
# 是否智能机                    -7.339148e-01 -2.863278e-01
# 合约客户                      -5.360059e+00 -3.608723e+00
# 合约类型                       4.118110e-01  1.763589e+00
# 账户余额                      -1.383972e-04 -6.987434e-05
# 月末欠费金额                  -7.024009e-02 -6.280929e-02
# 近三个月欠费次数              -2.314914e+00 -1.968475e+00
# 近三个月缴费金额               8.127753e-05  1.184311e-04
# 赠送费用                      -8.234064e-01 -7.986561e-01
# 赠送费用占出账收入比例        -5.302110e+00 -2.805590e+00
# 总计费时长                     3.515345e-03  5.938705e-03
# 被叫计费时长                  -7.456061e-03 -2.855837e-03
# 漫游通话次数                   1.016099e-03  5.629086e-03
# 主叫通话次数占比              -1.214744e+00 -4.622226e-01
# 上网计费流量                   2.637749e-04  3.658005e-04
# 出帐收入下降比例.比上月.       6.963364e-04  2.018112e-03
# 当月0通话天数                  4.925386e-02  7.837114e-02
# 连续3月平均出账收入            1.658800e-01  1.776232e-01
# 连续3月平均总通话次数          7.852762e-03  2.873983e-02
# 连续3月平均主叫次数           -2.733215e-02 -5.754663e-03
# 连续3月平均主叫计费时长       -7.481905e-03 -4.513359e-03
# 连续3月平均被叫次数           -3.069080e-02 -8.670937e-03
# 连续3月平均点对点短信发送条数 -8.935840e-03 -2.410540e-04
# 连续3月平均上网计费流量       -3.361233e-04 -2.348343e-04
# 近三月出账收入波动率           4.429877e+00  4.759986e+00
# 是否集团拍照                  -2.526148e+00 -1.201180e+00
# X12个月累计出账收入            3.080105e-03  3.854691e-03
# 当月是否停机                  -8.843691e-01 -4.191271e-01
# 本网通话时长占比               4.922301e-01  1.484347e+00
# 连续3月平均上网时长.分.       -6.070574e-05 -3.491796e-05
# 计费原始出账收入               8.057814e-01  8.148304e-01
# 合约剩余月数                   1.722746e-03  1.872472e-03
#-----------------------------
#使用前向逐步回归法再试试
min.model<-lm(出账收入~1,data=trainData)
fwdmodel<-step(min.model,direction="forward",scope = "~是否融合业务 + 用户当前积分余额 +
                   是否智能机 +  合约客户 +
               合约类型 + 账户余额 + 月末欠费金额 + 近三个月欠费次数 + 近三个月缴费金额 +
               赠送费用 + 赠送费用占出账收入比例 + 总计费时长 +
               被叫计费时长 + 漫游通话次数 + 主叫通话次数占比 + 上网计费流量 +
               出帐收入下降比例.比上月. + 当月0通话天数 + 连续3月平均出账收入 +
               连续3月平均总通话次数 + 连续3月平均主叫次数 + 连续3月平均主叫计费时长 +
               连续3月平均被叫次数 + 连续3月平均点对点短信发送条数 + 连续3月平均上网计费流量 +
               近三月出账收入波动率 +
               是否集团拍照 + X12个月累计出账收入 + 当月是否停机 + 本网通话时长占比 +
               连续3月平均上网时长.分. + 计费原始出账收入 + 合约剩余月数",data=trainData)
summary(fwdmodel)

#-----------------------------
#11、比较模型AIC
AIC(fwdmodel,bwdmodel)
#          df      AIC
# fwdmodel 32 281100.4
# bwdmodel 35 281093.0
#-----------------------------
#还可以利用以下函数查看变量重要性
relweights <-
  function(fit,...){                        
    R <- cor(fit$model)  
    nvar <- ncol(R)         
    rxx <- R[2:nvar, 2:nvar]
    rxy <- R[2:nvar, 1]     
    svd <- eigen(rxx)       
    evec <- svd$vectors                          
    ev <- svd$values        
    delta <- diag(sqrt(ev)) 
    lambda <- evec %*% delta %*% t(evec)       
    lambdasq <- lambda ^ 2  
    beta <- solve(lambda) %*% rxy          
    rsquare <- colSums(beta ^ 2)                  
    rawwgt <- lambdasq %*% beta ^ 2   
    import <- (rawwgt / rsquare) * 100
    lbls <- names(fit$model[2:nvar])  
    rownames(import) <- lbls
    colnames(import) <- "Weights"
    #import<-import[order(import["Weights"],decreasing=TRUE)]
    View(import)
    barplot(t(import),names.arg=lbls,
            ylab="% 重要性",
            xlab="预测变量",
            width=0.2,
            cex.axis = 0.8,
            cex.names = 0.6,
            las=2,
            main="预测变量的相对重要性",
            sub=paste("R-Square=", round(rsquare, digits=3)),
            ...) 
    return(import)
  }

relweights(bwdmodel, col="blue")

#-----------------------------
#12、使用模型做一下预测
options(digits = 3)
predict(bwdmodel, testData, interval = "prediction", level = 0.95)
# 得到预测结果与区间（fit预测值，lwr下限，upr上限）
#             fit       lwr       upr
# 1      1.04e+02  8.44e+01    124.54
# 7      9.73e+01  7.72e+01    117.45
# 8      4.22e+01  2.21e+01     62.24
# 13     1.31e+00 -1.88e+01     21.38
# 15     1.66e+01 -3.44e+00     36.68
# 17     9.28e+01  7.28e+01    112.91
# 18     1.89e+01 -1.13e+00     38.99
# 19     1.67e+01 -3.35e+00     36.77
# 20     4.01e+01  2.01e+01     60.20
#可以进一步将预测结果添加到数据框中，与实际值进行对比分析（自行练习）
#-----------------------------
#-----------------------------
#三、逻辑回归
#1、选择目标变量为是否离网，设置模式为二分类binomial
glmmodel<-glm(是否离网~.,data=trainData,family = binomial)
summary(glmmodel)
#2、根据统计量可以选择筛选相关指标重建模型
glmmodel2<-glm(是否离网~是否集团客户+是否融合业务+终端是否识别+入网渠道类型+靓号等级+信用额度
                   +月末欠费金额
                   +近三个月欠费次数
                   +赠送费用占出账收入比例
                   +主叫通话次数占比
                   +漫游通话次数占比
                   +交往圈个数
                   +最近一次通话距离天数
                   +近三天通话次数
                   +近三月通话次数波动率
                   +近三月主叫通话次数波动率
                   +近三月流量使用波动率
                   +近三月出账收入波动率
                   +X12个月累计出账收入
                   +当月是否停机
                   +本网通话时长占比
                   +连续3月平均上网时长.分.
                   +下行五位银行短信条数
                   +下行五位保险短信条数
                   +下行五位交通短信条数
                   +上行五位短号码条数
                   +识别出的短号码个数
                   +合约剩余月数,data=trainData,family = binomial)
#3、调用刚才的函数查看一下变量的重要性
relweights(glmmodel2,col="blue")

#可以根据重要性进一步筛选指标并重建逻辑回归模型
#-----------------------------
#4、分析查全率与查准率
#-----------------------------
#先看说明
#                       预测值
#                      0        1
# 实际值  0     a        b
#               1     c        d

# 其中，d是“实际为1而预测为1”的样本个数，c是“实际为1而预测为0”的样本个数，其余依此类推。
# 显然地，主对角线所占的比重越大，则预测效果越佳，这也是一个基本的评价指标——总体准确率(a+d)/(a+b+c+d)。
# 准确（分类）率=正确预测的正反例数/总数      Accuracy=(a+d)/(a+b+c+d)
# 误分类率=错误预测的正反例数/总数        Error rate=(b+c)/(a+b+c+d)=1-Accuracy
# 正例的覆盖率=正确预测到的正例数/实际正例总数
# Recall(True Positive Rate，or Sensitivity)=d/(c+d)
# 正例的命中率=正确预测到的正例数/预测正例总数
# Precision(Positive Predicted Value,PV+)=d/(b+d)
# 负例的命中率=正确预测到的负例个数/预测负例总数
# Negative predicted value(PV-)=a/(a+c)
# 5、ROC曲线
# 当预测效果较好时，ROC曲线凸向左上角的顶点。平移图中对角线，与ROC曲线相切，可以得到TPR较大而FPR较小的点。模型效果越好，则ROC曲线越远离对角线，极端的情形是ROC曲线经过（0，1）点，即将正例全部预测为正例而将负例全部预测为负例。ROC曲线下的面积可以定量地评价模型的效果，记作AUC，AUC越大则模型效果越好。
# 当我们分类的目标是将正例识别出来时，我们关注TPR，此时ROC曲线是评价模型效果的准绳。
#-----------------------------

#对模型做出预测结果
pre=predict(glmmodel2,type='response')
#将预测概率pre和实际结果放在一个数据框中
predata=data.frame(prob=pre,obs=trainData$是否离网)
#将预测概率按照从低到高排序
predata=predata[order(predata$prob),]
n=nrow(predata)
tpr=fpr=rep(0,n)
#根据不同的临界值threshold来计算TPR和FPR，之后绘制成图
for (i in 1:n){
 threshold=predata$prob[i]
 tp=sum(predata$prob>threshold&predata$obs==1)
 fp=sum(predata$prob>threshold&predata$obs==0)
 tn=sum(predata$prob)
 fn=sum(predata$prob)
 tpr[i]=tp/(tp+fn) #真正率
 fpr[i]=fp/(tn+fp) #假正率
}
plot(fpr,tpr,type='l')
abline(a=0,b=1)

#还可以使用专门绘制ROC曲线的包
install.packages("ROCR")
library(ROCR)
pred=prediction(pre,trainData$是否离网)
performance(pred,'auc')@y.values
perf=performance(pred,'tpr','fpr')
plot(perf)
#或者使用pROC包来绘制
install.packages("pROC")
library(pROC)
modelroc=roc(trainData$是否离网,pre)
plot(modelroc,print.auc=TRUE,auc.polygon=TRUE,grid=c(0.1,0.2),grid.col=c("green","red"),max.auc.polygon=TRUE,auc.polygon.col="blue",print.thres=TRUE)
#统计预测值与实际值
preldata<-within(predata,{
  prel<-NA
  prel[prob<0.5]<-0
  prel[prob>=0.5]<-1
})
table(preldata$obs,preldata$pre)
#                   0               1
# 0  16659(a)     1550(b)
# 1     1373(c)  17704(d)
# 查准率：d/(b+d)=17704/(17704+1550)=0.919497247
# 覆盖率：d/(c+d)=17704/(17704+1373)=0.928028516