广告投入是怎样提高新用户数的（岭回归及主成分回归） | R语言商业分析实践3

作者：杨奉山  R语言中文社区专栏作者，R语言小学生，和我一起为成为优秀的商业数据分析师努力吧 知乎专栏：https://zhuanlan.zhihu.com/YFSbda
最后一次重申：一个好的商业分析项目从来都不是从数据出发的，而是根据现象提出问题，之后根据问题从大局出发进行总体的（非数据上）分析（要清晰的把握自己的问题是什么，如何通过数据分析解决问题），根据这一步的分析建立具体数据分析框架，凭此再去确定要收集的数据和加工数据，最后分析得出结论解决问题。特别的是，在我们模型通过统计意义上的检验后，仍需结合实际情况推断模型的合理性；在模型未通过统计意义上的检验时，亦可结合实际找到出问题的地方。所以本章我们能看到如何结合实际对看似合格模型进行判断，同时结合实际情况找到问题的来源。（P.S.本节内容不涉及数据操作，主要包括多元回归及诊断、逐步回归、岭回归，主成分回归，多元回归的拟合三维图。可看只自己感兴趣的内容）第一部分：项目及思路分析

手游企业的广告部对投放的广告是否吸引了新玩家的进入产生了质疑，也想弄明白自己两大广告投放领域（电视及杂志）哪一个更有用。这一次的项目总体分析很简单，找出新玩家安装数目与广告投入额度的关系。分析框架有了，按这个框架我们初步决定收集的数据有：近10个月的电视与杂志上投入的广告费及期间的新用户数。

第二部分：数据分析

1.读取数据

(ad.data <- read.csv("ad_result.csv", header = T, stringsAsFactors = F))
    month  tvcm magazine install    
    1  2013-01  6358     5955   53948
    2  2013-02  8176     6069   57300
    3  2013-03  6853     5862   52057
    4  2013-04  5271     5247   44044
    5  2013-05  6473     6365   54063 6  2013-06  7682     6555   58097 7  2013-07  5666     5546   47407 8  2013-08  6659     6066   53333 9  2013-09  6066     5646   49918 10 2013-10 10090     6545   59963

2.画图观察install与其他两个变量是否存在相关性。

#观察有无相关性
ggplot(ad.data, aes(x = tvcm, y = install)) + geom_point()
ggplot(ad.data, aes(x = magazine, y = install)) + geom_point()

从图上看存在很强的正相关性。3.0初步的回归分析及结合现实的统计诊断（本章重点1）

#初步回归
fit <- lm(install ~ ., data = ad.data[, c("install", "tvcm", "magazine")])
summary(fit)

我们对回归结果做简单分析：总体F检验的p值（5.967e-05）很小，但是截距项及tvcm的p值很大，特别是截距项p值已经高达0.98，说明截距项极不显著（在计量里直接删除截距项做回归即可）。但我们现在应该结合实际分析：截距项相当于与五位数的install而言可以忽略不计，那么说明了什么？说明install与另外两个变量的回归是一个原点的曲面，这在现实里是不科学的（否则公司只要把广告宣传做好就啥都不用干了），因此可以初步推测两个变量tvcm和magazine存在一定程度相关性（被过度重视）。

(R <- cor(ad.data[,2:3]))#相关系数判断
#0.7723106。存在一定相关性

3.1逐步回归法（变量过少不适用本方法，因此上来笔者就失败了，但是数据分析过程本就不会一帆风顺）

step(fit)

#结果当然是毫无改变
Start: AIC=147.14
install ~ tvcm + magazine
Df Sum of Sq RSS AIC
<none> 13473540 147.14
- tvcm 1 13315165 26788705 152.01
- magazine 1 35310619 48784159 158.00
Call:
lm(formula = install ~ tvcm + magazine, data = ad.data[, c("install", 
"tvcm", "magazine")])
Coefficients:
(Intercept) tvcm magazine 
188.174 1.361 7.2503.2岭回归法：首先交叉验证选lamda值

#岭回归
fit.ridge_lambda <- cv.glmnet(alpha = 0,y=ad.data$install, x = as.matrix(ad.data[, c( "tvcm", "magazine")]))
fit.ridge_lambda$lambda.min

[1] 1213.469然后按lamda=1213.469进行岭回归：

fit.ridge <- glmnet(alpha = 0,y=ad.data$install, x = as.matrix(ad.data[, c( "tvcm", "magazine")]),lambda = 1213.469)
coef(fit.ridge)

得到系数参数值：
3 x 1 sparse Matrix of class "dgCMatrix"
s0
(Intercept) 9152.276952
tvcm 1.360070
magazine 5.753183对比之前的简单回归，截距扩大了50倍，也更能解释现实意义。3.3主成分回归（本章重点方法）

#第一步 主成分提取
pra <- prcomp(~tvcm+magazine,data=ad.data,scale=T)
summary(pra)
# Importance of components:
#   PC1    PC2
# Standard deviation     1.3313 0.4772
# Proportion of Variance 0.8862 0.1138
# Cumulative Proportion  0.8862 1.0000
ad.data.pra <- transform(ad.data,z1=pra$x[,1])

注意这里我们舍弃第二个成份，只用了第一个主成分。

#第一次对主成分回归
pra.lm0 <- lm(install ~ z1,data=ad.data.pra)
summary(pra.lm0)

这里看到，各个p值均很小，并且Adjusted R-squared也不错，但是我们还应对残差进行回归诊断。

plot(pra.lm0)

从残差图来看，我们还应该加上二次项来第二次回归。

# 第二次对主成分回归
pra.lm <- lm(install ~ z1+I(z1^2),data=ad.data.pra)
summary(pra.lm)

整体来看，比第一次主成分回归有一定的进步，下面我们再看残差图：

plot(pra.lm)

很正常到这里数据分析的工作就结束了，但是别忘了你的上司要的回归方程并不是关于主成分的，而是关于原变量的，因此还应该进行转换。笔者没有找到R语言中自动代回原变量的函数（如果有大佬知道请留言，不胜感激），因此用了笨的方法：手工推导，以第一次主成分回归为例（第二次主成分回归有二次项，笔者又不会用R语言对含未知数表达式的二次项开方，有知道如何对含未知数的二次项表达式开方也请留言）：

#转化系数
bet0 <- coef(pra.lm0)[1]-(sum(pra$rotation[1,1]*pra$center/pra$scale))*coef(pra.lm0)[2]
beta <- coef(pra.lm0)[2]*(1/pra$scale)*pra$rotation[1,1]
c(bet0,beta)
#(Intercept)        tvcm    magazine
#5598.356864    1.788569    5.850864

4.番外篇——三元回归的立体图这里笔者只给出第一次简单回归画图的代码（有兴趣的可以看《R语言可视化手册》：

#设置函数
predictgrid <- function(model, xvar, yvar, zvar, res = 16, type = NULL) {
 # Find the range of the predictor variable. This works for lm and glm
 # and some others, but may require customization for others.
 xrange <- range(model$model[[xvar]])
 yrange <- range(model$model[[yvar]])
 
 newdata <- expand.grid(x = seq(xrange[1], xrange[2], length.out = res),
                        y = seq(yrange[1], yrange[2], length.out = res))
 names(newdata) <- c(xvar, yvar)
 newdata[[zvar]] <- predict(model, newdata = newdata, type = type)
 newdata
}

df2mat <- function(p, xvar = NULL, yvar = NULL, zvar = NULL) {
 if (is.null(xvar)) xvar <- names(p)[1]
 if (is.null(yvar)) yvar <- names(p)[2]
 if (is.null(zvar)) zvar <- names(p)[3]
 
 x <- unique(p[[xvar]])
 y <- unique(p[[yvar]])
 z <- matrix(p[[zvar]], nrow = length(y), ncol = length(x))
 
 m <- list(x, y, z)
 names(m) <- c(xvar, yvar, zvar)
 m
}

# Function to interleave the elements of two vectors
interleave <- function(v1, v2)  as.vector(rbind(v1,v2))
###设置参数
ad.data$pred <- predict(fit)
install_df<- predictgrid(fit,'tvcm','magazine','install')
install_list <- df2mat(install_df)

#为了一次性成图，最好加上{}
{
 plot3d(ad.data$tvcm,ad.data$magazine,ad.data$install,
      size=0.5,type='s',axes=F,
      xlab = '',ylab='',zlab='')

spheres3d(ad.data$tvcm,ad.data$magazine,ad.data$pred,type='s',size=0.5)

segments3d(interleave(ad.data$tvcm,ad.data$tvcm),interleave(ad.data$magazine,ad.data$magazine),
      interleave(ad.data$install,ad.data$pred),col='blue',alpha=0.4)

surface3d(install_list$tvcm,install_list$magazine,install_list$install,alpha=0.4,front='lines',back='lines')

rgl.bbox(color='grey50',emission='grey50',xlen=0,ylen=0,zlen=0)

rgl.material(color='black')

axes3d(edges=c("x--", "y+-", "z--"),
      ntick=6,                       # Attempt 6 tick marks on each side
      cex=.75)                       # Smaller font

# Add axis labels. 'line' specifies how far to set the label from the axis.
mtext3d("TV",edge="x--", line=2)
mtext3d("Magazine",edge="y+-", line=3)
mtext3d("Install",edge="z--", line=3)
}
#动态观察
play3d(spin3d())

公众号后台回复关键字即可学习回复 R               R语言快速入门免费视频 
回复 统计          统计方法及其在R中的实现
回复 用户画像   民生银行客户画像搭建与应用 
回复 大数据      大数据系列免费视频教程
回复 可视化      利用R语言做数据可视化
回复 数据挖掘   数据挖掘算法原理解释与应用
回复 机器学习   R&Python机器学习入门