2.3 R商業分析3：廣告投入是怎樣提高新用戶數的（嶺回歸及主成分回歸）

05-13

（本文除數據及番外篇以外所有內容包括代碼均為原創，轉載請註明出處：知乎專欄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.250

3.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.0000ad.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 vectorsinterleave <- 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語言中對含未知數的二次項表達式開方（如果有大佬知道也請留言，不勝感激）