R在市场调查中的应用--逻辑回归和LDA

本文链接：https://blog.csdn.net/JohnSuperD/article/details/77484670

一、分类模型的应用场景

我们一般接触的统计学习的模型多半是因变量Y为数值型的模型，这一类问题可以统称为回归（regression）问题，当Y不再是数值而是类别的时候，就需要使用类似逻辑回归、LDA之类的统计模型，而这一类的问题我们统称为分类（classification）问题。
逻辑回归一般应用于Y只有两个值的情况，比如Yes或者No。在市场调查的应用场景下，Y可以是再次购买与否。而LDA则可应用于Y有2个值及以上的情形。
不同于常见的线性回归模型依赖于最小二乘法来拟合最适模型，逻辑回归依靠最大熵（Maxium entropy）来确定。可以这样地理解用最大熵来拟合模型的思想–即最大熵模型就是在给定的条件下，概率最高、可能性最大的模型。具体可以参见《Statistical Rethinking》第九章。
LDA（Linear Discriminant Analysis）恰如其名，其实是一种线性判定的方式，即该模型给定的判定边界是线性的。这里有一个关于深入理解LDA算法的教程，
下面我们将结合实际的市场调查中的应用来理解逻辑回归和LDA的应用。

二、数据的准备

现在我们来看探讨一个影响顾客购买主题公园年卡行为的案例–研究促销的渠道（包括通过邮件、纸质信件、在公园的促销）和是否采用捆包销售对客户购买年卡的行为是否会产生影响。这里的数据来自《R for Marketing Research and Analytics》第九章。

library(readr)
data_pass<- read_csv("http://r-marketing.r-forge.r-project.org/data/rintro-chapter9.csv")

str(data_pass)

可以看到，data_pass记录了客户购买主题公园年卡（pass）的情况，pass变量在这里就是我们响应变量，是包含了“YesPass”和“Nopass”的二值变量。channel变量记录了客户的购买年卡的渠道，包含了电子邮件、信封，公园三种方式。而bundle变量则记录了是不是采取了捆包销售的方式，比如和“免费停车”服务捆绑起来。
首先我们可以用table()函数来看一下数据的分布状况。

table(data_pass)

输出结果：

, , Pass = NoPass

       Promo
Channel Bundle NoBundle
  Email     83      485
  Mail     449      278
  Park     223       49

, , Pass = YesPass

       Promo
Channel Bundle NoBundle
  Email     38       27
  Mail     242      359
  Park     639      284

因为有三个维度，而且作为Y变量对Pass变量被安排在两个表中，所以很难直观地看出各因素对Pass的影响。因此可以使用vcd包的doubledecker函数来进行更为直观的呈现。

library(vcd)
doubledecker(table(data_pass))

这里写图片描述
如何看这张图呢？比如在相同的channel的情况下，可以看到除了在使用Email渠道以外，捆绑销售的效果不如没有捆包的效果好，直接体现在“Nobundle”的“YesPass”的黑色部分的比例要大于“Bundle”来的大。但是，是不是就能断言Bundle这种形式没有任何的促销效果呢？还不能就此断言，以Park渠道为例，可以看到Bundle和NoBundle的数量相差极大，所以光看“YesPass”和“NoPass”的比例是不可以的。

三、逻辑回归（logistic regression）的应用

data_pass$Promo<-factor(data_pass$Promo,levels=c("NoBundle","Bundle"))
#改变Promo变量中level的顺序，这样“捆绑销售”将作为dummy变量，可以直接从后面模型的系数中看到它与预测变量之间的关系。

pass_fact<-as.data.frame(lapply(data_pass,as.factor))
#将所有地变量都变味factor

使用glm()来拟合逻辑回归模型，注意将其中的family参数设置为binomial

pass.LR1<-glm(Pass~ Promo+Channel,data = pass_fact,family = binomial())
summary(pass.LR1)

输出结果：

Call:
glm(formula = Pass ~ Promo + Channel, family = binomial(), data = pass_fact)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9079  -0.9883   0.5946   0.7637   2.3272  

Coefficients:
            Estimate Std. Error z value
(Intercept) -2.07860    0.13167 -15.787
PromoBundle -0.56022    0.09031  -6.203
ChannelMail  2.17617    0.14651  14.854
ChannelPark  3.72176    0.15964  23.313
            Pr(>|z|)    
(Intercept)  < 2e-16 ***
PromoBundle 5.54e-10 ***
ChannelMail  < 2e-16 ***
ChannelPark  < 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: 4375.0  on 3155  degrees of freedom
Residual deviance: 3490.2  on 3152  degrees of freedom
AIC: 3498.2

Number of Fisher Scoring iterations: 4

结合前面的模型概况中的coefficients一项以及每一项的Pr(>|z|)，我们可以得到不采用捆绑销售的方法似乎反而更能促进年卡的销售，另外相较于电子邮件，纸质邮件和在公园的促销活动更有效果，特别是公园的促销活动，非常有效果。
但是是不是捆绑式的促销活动真没有效果呢？从上面的模型来看似乎是这样，但是我们必须要考虑到的是预测变量之间的相关性，以及预测变量间的交互作用对响应变量地影响。
首先我们可以通过开方检验来查看两个预测变量之间是不是存在着某种相关性。

tbl1<-table(pass_fact$Channel,pass_fact$Promo)
tbl1

chisq.test(tbl1)

输出结果：

        NoBundle Bundle
  Email      512    121
  Mail       637    691
  Park       333    862

    Pearson's Chi-squared test

data:  tbl1
X-squared = 467.94, df = 2, p-value < 2.2e-16

没错，从p-value看，Promo变量和Channel变量存在着相关性（主要体现在respondets的数量上，可以看到Park&Bundle的数量明显非常多），而这种相关性会不会影响到最后的模型呢？
为了一探究竟，我们可以先考虑一个只有Promo一个变量的逻辑回归模型。

pass.LR2<-glm(Pass~ Promo,data = pass_fact,family = binomial())
summary(pass.LR2)

输出结果：

Call:
glm(formula = Pass ~ Promo, family = binomial(), data = pass_fact)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.262  -1.097   1.095   1.095   1.260  

Coefficients:
            Estimate Std. Error z value
(Intercept) -0.19222    0.05219  -3.683
PromoBundle  0.38879    0.07167   5.425
            Pr(>|z|)    
(Intercept) 0.000231 ***
PromoBundle 5.81e-08 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 4375.0  on 3155  degrees of freedom
Residual deviance: 4345.4  on 3154  degrees of freedom
AIC: 4349.4

Number of Fisher Scoring iterations: 3

可以看到PromoBundle的系数是正的，而且从p值上来看也是有意的。这样，我们得到了和之前的多变量模型完全不一样的结果，为甚么会这样？结合前面的开放检验的结果就说的通了，所以当计算Park对响应变量产生影响的时候，其实也包含了Bundle的效果。因此，有必要在模型中引入交互项。

pass.LR3<-glm(Pass~ Promo*Channel,data = pass_fact,family = binomial)

summary(pass.LR3)

输出结果：

Call:
glm(formula = Pass ~ Promo * Channel, family = binomial, data = pass_fact)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9577  -0.9286   0.5642   0.7738   2.4259  

Coefficients:
                        Estimate Std. Error
(Intercept)              -2.8883     0.1977
PromoBundle               2.1071     0.2783
ChannelMail               3.1440     0.2133
ChannelPark               4.6455     0.2510
PromoBundle:ChannelMail  -2.9808     0.3003
PromoBundle:ChannelPark  -2.8115     0.3278
                        z value Pr(>|z|)    
(Intercept)             -14.608  < 2e-16 ***
PromoBundle               7.571 3.71e-14 ***
ChannelMail              14.743  < 2e-16 ***
ChannelPark              18.504  < 2e-16 ***
PromoBundle:ChannelMail  -9.925  < 2e-16 ***
PromoBundle:ChannelPark  -8.577  < 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: 4375.0  on 3155  degrees of freedom
Residual deviance: 3393.5  on 3150  degrees of freedom
AIC: 3405.5

Number of Fisher Scoring iterations: 5

这样捆绑销售和在公园的促销活动都被认为是对促进年卡销售的活动。只不过这两者的交互效果却是负的，也就是在公园销售是同时采用捆绑销售的做法是不利于销售的（因为这只是假想数据未必符合现实）。
最后，我们还可以使用anova以及AIC函数来比较不同的模型的拟合效果。

anova(pass.LR1,pass.LR3)

AIC(pass.LR1,pass.LR3)

输出结果：

Analysis of Deviance Table

Model 1: Pass ~ Promo + Channel
Model 2: Pass ~ Promo * Channel
  Resid. Df Resid. Dev Df Deviance
1      3152     3490.2            
2      3150     3393.5  2   96.661

#AIC

df
<dbl>
AIC
<dbl>
pass.LR1    4   3498.209        
pass.LR3    6   3405.547

可以看到引入了交互项的模型，有着更低的deciance（少了96.661），同时AIC值也更小，所以可以把pass.LR3看作为更好的模型。

为了验证模型的拟合效果，往往需要将模型应用于新的数据。当然在没有新数据的情况下可以先用将原数据一分为二的办法，即随机地将一半数据分配为训练子集，用以训练模型，另一部分则用来验证模型的拟合性。这样做的好处是可以方便多个模型之间的比较。

set.seed(020)

train<-sample(1:nrow(pass_fact),0.5*nrow(pass_fact))

pass.LR4<-glm(Pass~ Promo*Channel,data = pass_fact[train,],family = binomial)#

把数据随机地分为train set和test set（即非train）之后，我们用这两个数据集分别地拟合以及验证模型的拟合效果。

LR_props<-predict(pass.LR4,data=pass_fact[-train,],type = "response") #注意这里的type参数需要设置为“response”

LR_pred<-rep('NoPass',length(LR_props))
LR_pred[LR_props>0.5]<-"YesPass"

#correct rate
mean(pass_fact[-train,]$Pass==LR_pred)
#[1] 0.486692

prop.table(table(LR_pred,pass_fact[-train,]$Pass),2)
#Ture positive rate is :0.562

可以看到，该模型整体的预测效果并不是很理想，只预测对了48.7%，和瞎蒙差不多了，只不过TRUE positive rate还是达到了56.2%，也就是说该模型成功预测了56%购买了年卡的消费者的情况。

四、线性判别方法（LDA）的应用

接下来我们看一下LDA方法的拟合效果。为了方便比较，我们也使用训练集来拟合模型，用验证集来确认模型的拟合效果。

library(MASS)
pass.LDA1<-lda(Pass~Channel*Promo,data=pass_fact[train,])
pass.LDA1

输出结果：

Call:
lda(Pass ~ Channel * Promo, data = pass_fact[train, ])

Prior probabilities of groups:
  NoPass  YesPass 
0.495564 0.504436 

Group means:
        ChannelMail ChannelPark PromoBundle
NoPass    0.4539642   0.1828645   0.4833760
YesPass   0.3756281   0.5791457   0.5979899
        ChannelMail:PromoBundle
NoPass                0.2813299
YesPass               0.1620603
        ChannelPark:PromoBundle
NoPass                0.1457801
YesPass               0.4120603

Coefficients of linear discriminants:
                              LD1
ChannelMail              2.254463
ChannelPark              3.464286
PromoBundle              1.080533
ChannelMail:PromoBundle -1.942197
ChannelPark:PromoBundle -1.442668

接下来用该模型进行对test set的预测

LDA_pre<-predict(pass.LDA1,data=pass_fact[-train,])
LDA_class<-LDA_pre$class

mean(LDA_class==pass_fact[-train,]$Pass)
# 0.486692

prop.table(table(LDA_class,pass_fact[-train,]$Pass),2)

#True positive
0.5624

不要怀疑你的眼睛，不论是逻辑回归还是LDA，这两者的预测结果是一模一样的！
为甚会这样？
事实上，不论是逻辑回归还是LDA都是线性判别方法，它们的Decision Boundary（判别边界）是很相似的，通过比较pass.LR4和pass.LDA1的系数即可。另外一点就是数据本身的问题，因为pass_fact的预测变量都是factor，而且level也不多（只有2*3六种可能性），所以大多数的数据其实都是重叠在一起的，如果预测变量中又一个数值型数据的话，结果或许会大不相同。另外，因为预测变量是factor类型，所以也不能使用kmeans(需要计算数据点之间的距离)。