# 3.3.2 逻辑回归

    ##### 3.3 回归分析 #####
    rm(list = ls())  # 清空工作空间

    #### 3.3.2 逻辑回归 ####
    # 加载所需R包
    # install.packages("ggplot2")
    library(ggplot2)
    # 读入数据
    dat1 = read.csv("JuneTrain.csv")
    dat2 = read.csv("JulyTest.csv")
    head(dat1)  # 查看数据的前几行
    ##          sbid tenure expense COUNT perperson  entropy   chgexpense
    ##  1 1445479416   3608  100.39    48  6.844792 2.931435 -0.003375360
    ##  2 1445479496   3596   50.54    38  8.485000 2.075424 -0.132211538
    ##  3 1445479746   3604   29.49    85  4.537882 3.702529 -0.829438982
    ##  4 1445480426   2252   97.29    97  7.127010 2.773227  0.025616698
    ##  5 1445480735   2155  228.77    82 11.652439 3.096718  0.044611872
    ##  6 1445480785   2158   95.25   144  2.998958 4.507604 -0.002095338
    ##       chgcount churn
    ##  1 -0.21311475     0
    ##  2 -0.32142857     0
    ##  3 -0.27350427     0
    ##  4  0.25974026     0
    ##  5  0.03797468     0
    ##  6  0.03597122     0
    # 重新命名数据以及排列
    dat1 = dat1[, c("tenure", "expense", "COUNT", "perperson", "entropy", "chgexpense", "chgcount", "churn")]
    colnames(dat1) = c("tenure", "expense", "count", "perperson", "entropy", "chgexpense", "chgcount", "churn")
    dat2 = dat2[, c("tenure", "expense", "COUNT", "perperson", "entropy", "chgexpense", "chgcount", "churn")]
    colnames(dat2) = c("tenure", "expense", "count", "perperson", "entropy", "chgexpense", "chgcount", "churn")

    # churn:是否流失;tenure:在网时长;expense:当月话费;count:通话人数;
    # perperson:人均通话时长;entropy:通话时长分布;chgexpense:花费变化率;chgcount:通话人数变化率

    # 描述分析，用ggplot2画对比箱线图。
    # 描述分析内容 "tenure"     "expense"    "count"      "perperson"  "entropy"    "chgexpense" "chgcount"
    # tenure和是否流失
    p1 = ggplot(dat1, aes(x = as.factor(churn), y = tenure, fill = as.factor(churn)))  + geom_boxplot() +
      guides(fill = FALSE) + theme_minimal() + xlab("是否流失") + ylab("在网时长") +
      theme(axis.title.x = element_text(size = 16, face = "bold"),
            axis.text.x = element_text(size = 12, face = "bold"),
            axis.title.y = element_text(size = 16, face = "bold"),
            axis.text.y = element_text(size = 12, face = "bold")) +
      scale_fill_hue(c = 45, l = 80)
    # expense和是否流失
    p2 = ggplot(dat1, aes(x = as.factor(churn), y = expense, fill = as.factor(churn)))  + geom_boxplot() +
      guides(fill = FALSE) + theme_minimal() + xlab("是否流失") + ylab("当月话费") +
      theme(axis.title.x = element_text(size = 16, face = "bold"),
            axis.text.x = element_text(size =12, face = "bold"),
            axis.title.y = element_text(size = 16, face = "bold"),
            axis.text.y = element_text(size = 12, face = "bold")) +
      scale_fill_hue(c = 45, l = 80)
    # 加载所需R包
    # install.packages("gridExtra")
    library(gridExtra)
    # 将2张图并列输出
    grid.arrange(p1,p2,ncol = 2)

    # count和是否流失
    p3 = ggplot(dat1, aes(x = as.factor(churn), y = count, fill = as.factor(churn)))  + geom_boxplot() +
      guides(fill = FALSE) + theme_minimal() + xlab("是否流失") + ylab("通话人数") +
      theme(axis.title.x = element_text(size = 16, face = "bold"),
            axis.text.x = element_text(size = 12, face = "bold"),
            axis.title.y = element_text(size = 16, face = "bold"),
            axis.text.y = element_text(size = 12, face = "bold")) +
      scale_fill_hue(c = 45, l = 80)
    # perperson和是否流失
    p4 = ggplot(dat1, aes(x = as.factor(churn), y = perperson, fill = as.factor(churn))) + geom_boxplot() +
      guides(fill = FALSE) + theme_minimal() + xlab("是否流失") + ylab("人均通话时长") +
      theme(axis.title.x = element_text(size = 16, face = "bold"),
            axis.text.x = element_text(size = 12, face = "bold"),
            axis.title.y = element_text(size = 16, face = "bold"),
            axis.text.y = element_text(size = 12, face = "bold")) +
      scale_fill_hue(c = 45, l = 80)
    # 将2张图并列输出
    grid.arrange(p3, p4, ncol = 2)

    # entropy和是否流失
    p5 = ggplot(dat1, aes(x = as.factor(churn), y = entropy, fill = as.factor(churn)))  + geom_boxplot() +
      guides(fill = FALSE) + theme_minimal() + xlab("是否流失") + ylab("通话时长分布") +
      theme(axis.title.x = element_text(size = 16, face = "bold"),
            axis.text.x = element_text(size = 12, face = "bold"),
            axis.title.y = element_text(size = 16, face = "bold"),
            axis.text.y = element_text(size = 12, face = "bold")) +
      scale_fill_hue(c = 45, l = 80)
    # chgcount和是否流失
    p6 = ggplot(dat1, aes(x = as.factor(churn), y = chgcount, fill = as.factor(churn))) + geom_boxplot() +
      guides(fill = FALSE) + theme_minimal() + xlab("是否流失") + ylab("通话人数变化率") +
      theme(axis.title.x = element_text(size = 16, face = "bold"),
            axis.text.x = element_text(size = 12, face = "bold"),
            axis.title.y = element_text(size = 16, face = "bold"),
            axis.text.y = element_text(size = 12, face = "bold")) +
      scale_fill_hue(c = 45, l = 80)
    # 将2张图并列输出
    grid.arrange(p5, p6, ncol = 2)

    # chgexpense和是否流失
    p7 = ggplot(dat1, aes(x = as.factor(churn), y = chgexpense, fill = as.factor(churn))) + geom_boxplot() +
      guides(fill = FALSE) + theme_minimal() + xlab("是否流失") + ylab("花费变化率") +
      theme(axis.title.x = element_text(size = 16, face ="bold"),
            axis.text.x = element_text(size = 12, face = "bold"),
            axis.title.y = element_text(size = 16, face = "bold"),
            axis.text.y = element_text(size = 12, face = "bold")) +
      scale_fill_hue(c = 45, l = 80)

    # 计算异常值
    flag = apply(dat1[, 1:(ncol(dat1) - 1)], 2, function(s) {
    # s = log(s - min(s) + 1)
      m1 = mean(s)
      s1 = sd(s)
      res = (s > m1 - 3 * s1 & s < m1 + 3 * s1) * 1
      return(res)
    })
    flagsub = apply(flag, 1, function(s) {mean(s) == 1})
    # 取出测试集中无异常值的部分
    dat1 = dat1[flagsub, ]

    # 标准化，用训练集的标准标准化测试集
    for(i in 1:(ncol(dat1) - 1))
    {
      mm1 = mean(dat1[, i])
      ss1 = sd(dat1[, i])
      dat1[, i] = (dat1[, i] - mm1) / ss1
      dat2[, i] = (dat2[, i] - mm1) / ss1
    }

    # 建立模型
    lm1 = glm(churn ~., data = dat1, family = binomial())
    summary(lm1)
    ##
    ##  Call:
    ##  glm(formula = churn ~ ., family = binomial(), data = dat1)
    ##
    ##  Deviance Residuals:
    ##      Min       1Q   Median       3Q      Max
    ##  -0.9070  -0.1906  -0.1357  -0.0921   4.0309
    ##
    ##  Coefficients:
    ##              Estimate Std. Error z value Pr(>|z|)
    ##  (Intercept) -4.71208    0.06070 -77.624  < 2e-16 ***
    ##  tenure      -0.30973    0.05732  -5.404 6.52e-08 ***
    ##  expense     -0.19255    0.05246  -3.671 0.000242 ***
    ##  count       -0.52674    0.11092  -4.749 2.04e-06 ***
    ##  perperson   -0.14177    0.03997  -3.547 0.000390 ***
    ##  entropy     -0.26415    0.07356  -3.591 0.000330 ***
    ##  chgexpense  -0.19557    0.04304  -4.544 5.52e-06 ***
    ##  chgcount    -0.38860    0.04270  -9.101  < 2e-16 ***
    ##  ---
    ##  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ##  (Dispersion parameter for binomial family taken to be 1)
    ##
    ##      Null deviance: 6619.6  on 42623  degrees of freedom
    ##  Residual deviance: 5977.7  on 42616  degrees of freedom
    ##  AIC: 5993.7
    ##
    ##  Number of Fisher Scoring iterations: 8

    # 模型预测
    Yhat = predict(lm1, newdata = dat2, type = "response")
    ypre1 = 1 * (Yhat > 0.5)
    table(ypre1, dat2$churn)
    ##
    ##  ypre1     0     1
    ##      0 46001   447
    ypre2 = 1 * (Yhat > mean(dat2$churn))
    table(ypre2, dat2$churn)
    ##
    ##  ypre2     0     1
    ##      0 27242   108
    ##      1 18759   339

    # 覆盖率捕获率曲线
    sub = seq(0, 1, 1 / 7)
    tol = sum(dat2$churn)
    catch = sapply(sub, function(s) {
      ss = quantile(Yhat, 1 - s)
      res = sum(dat2$churn[Yhat > ss]) / tol
      return(res)
    })
    plot(sub, catch, type = "l", xlab = "覆盖率", ylab = "捕获率")

    # 系数图示
    coef = lm1$coefficients[-1]
    coef = sort(coef)
    barplot(coef, col = rainbow(10), width = 1)

    # 加载所需R包
    # install.packages("pROC")
    library(pROC)
    # 生成ROC曲线
    plot.roc(dat2$churn, Yhat, col = "red", lwd = 2, xaxs = "i", yaxs = "i")

    # 比较ROC曲线优劣
    lm2 = glm(churn ~ chgcount, data = dat1, family = binomial())
    Yhat2 = predict(lm2, newdata = dat2, type = "response")
    plot.roc(dat2$churn, Yhat2, col = "blue", lwd = 2, xaxs = "i", yaxs = "i")
    lines.roc(dat2$churn, Yhat, col = "red", lwd = 2)

    # auc曲线
    auc(dat2$churn, Yhat)
    ##  Area under the curve: 0.7554