用多元线性回归预测网页访问量(R语言)

前言

该问题来源于《机器学习:实用案例解析》中的第5章。在书中，已经对该问题给出了一种解决方案，但是我觉得写的还是太简略了一些，没有把考虑问题的整个思路给写出来，所以，在这里给出我的一些想法。

问题简述

我们的任务就是根据给定的数据集（TOP1000的互联网站数据）建立一个回归模型，然后根据任意给定的一组网站数据，预测出该网站的网页访问量。

解决方案

这里我们针对的是多元线性回归这个方法，并不是针对预测网页访问量这件事（只是拿这个数据集来练习一下~），看看多元线性回归能拟合到怎样的程度，所以如果结果不尽如人意也就不要纠结了~
好了，首先拿到数据集之后自然要看一下数据集是长什么样的。

topSites <- read.csv('top_1000_sites.tsv', sep='\t', stringsAsFactors = FALSE)
head(topSites)

得到

同样由于对“法”不对“事”，我们这里只考虑UniqueVisitors（用户数），HasAdvertising（网页是否有广告）和InEnglish（是否是英文网站）三者对PageViews（网页访问量）的影响。
先来看一下UniqueVisitors和PageViews大概是个什么关系

plot(topSites$PageViews, topSites$UniqueVisitors)

得到

图1：用户数和访问量的关系图
可以看出这并没有很好的线性关系，当然我们也可以强行来做一下线性拟合看下效果~

lm.sol <- lm(PageViews~UniqueVisitors, data = topSites)
lines(topSites$PageViews, fitted(lm.sol))

图2：有拟合线的用户数和访问量的关系
再来看下具体的结果

summary(lm.sol)

在这里我们主要看三个指标（详细的可参见参考文献[2]中第6章），一个是Coefficients后面的 * 的个数，越多说明显著性越高，当然也就越好；再看Residual standard error（残差的标准差），越小越好，这里大的有些吓人，这主要也和数据本身的数量级有些关系，但这里我们显然是不满意的；最后看Mutiple R-squared（相关系数的平方），取值为0~1之间，越靠近1越好，这里的0.3421只能说是有点线性相关，相关性不强。
好了，接下来我们来改进模型。
首先，通过观察数据不难发现，UniqueVisitors(x)和PageViews(y)有点幂函数的意思，说的更具体些，有点介于 $y=a \cdot x$ 和 $y=a \cdot x^2$ 之间的意思~
我们来试下看~
假设UniqueVisitors(x)和PageViews(y)有如下关系：
$y = a \cdot x^b$
对两边分别取对数有
$lny = lna + blnx$
这样一来 $lny$和 $lnx$就成了一个线性的关系了，然后就可以名称言顺地放进去线性拟合了，其中 $lna$为常数项， $b$为斜率。

plot(log(topSites$UniqueVisitors), log(topSites$PageViews))
lm.solLog <- lm(log(PageViews)~log(UniqueVisitors), data = topSites)
abline(lm.solLog)

图3：取对数后的线性拟合图
再来用数据说话看下~

summary(lm.solLog)

再来看下这三个值，显著性依旧比较理想，残差的标准差1.084也可以接受（主要是因为取对数之后数据的值的数量级都变小了），相关系数的平方0.4616也挺不错了，所以我们认为这个结果就线性回归而言已经比较理想了。
注意到这里估计出来的截距为-2.83441，斜率为1.33628，也就是说估计出来的方程为
$lny = -2.83441 + 1.33628 \cdot lnx$
即
$y = e^{-2.83441} \cdot x^{1.33628}$
看来我们的猜测还是比较准确的~
当然也可以猜测UniqueVisitors和PageViews是指数或者其他的关系，我都试了一下，不太理想，所以我们还是用幂函数的关系~
接下来就是另外两个变量HasAdvertising，InEnglish与PageViews的关系了。但遗憾的是这两个变量的值是YES和NO，这怎么线性拟合？。。。
这里涉及到一个虚拟变量（离散变量）的问题（详情可参见参考文献[1]第5章），我们来增加两列数据isEnglishYes和hasAdYes。

topSites[topSites$InEnglish == "Yes", "isEnglishYes"] <- 1
topSites[is.na(topSites$isEnglishYes), "isEnglishYes"] <- 0

topSites[topSites$HasAdvertising == "Yes", "hasAdYes"] <- 1
topSites[is.na(topSites$hasAdYes), "hasAdYes"] <- 0

这里把 InEnglish == “Yes” 和 HasAdvertising == “Yes” 的行置为1，InEnglish == “No” 和 HasAdvertising == “No” 的行置为0。这样做之后就可以对最终拟合结果的截距和斜率产生影响。
不用搞个 isEnglishNo 这样的列出来，因为 isEnglishNo 这个列是由 isEnglishYes 这列决定的，两个都放进去就会有一个多重共线性的问题。
比如下面这样就是对截距产生影响：

$lny = a + b \cdot lnx + c \cdot IsEnglish + d \cdot HasAdvertising$
也就是说是不是英语网站和有没有广告投入这几种情况下的截距是不同的，但是斜率是一样的。
又比如下面这样是对斜率产生影响：
$lny = a + b \cdot lnx + c \cdot isEnglish \cdot lnx + d \cdot HasAdvertising \cdot lnx$
也就是说是不是英语网站和有没有广告投入这几种情况下的斜率是不同的，但截距是一样的。
这就是虚拟变量的处理方法，我们下面来综合考虑一下上面两种情况：
$lny = a + b \cdot lnx + c \cdot IsEnglish + d \cdot HasAdvertising + e \cdot isEnglish \cdot lnx + f \cdot HasAdvertising \cdot lnx$

lm.solMul2 <- lm(log(PageViews)~log(UniqueVisitors)+hasAdYes+isEnglishYes+hasAdYes*log(UniqueVisitors)+isEnglishYes*log(UniqueVisitors), data = topSites)
summary(lm.solMul2)

可以得到

我的天！竟然只有一个参数估计的显著性是三颗 * 的，其它的连 . 都没有！也就是说没有通过参数检验，这个模型肯定用不了了，后面的看都不用看了。
没关系，这其实也是在意料之中的事，也是多元回归中经常碰到的问题。接下来，我们就用逐步回归来解决这个问题。
逐步回归分为：

• 向前引入法（forward）：逐渐增加变量，直到AIC值最优。
• 向后剔除法（backward）：逐渐删减变量，直到AIC值最优。
• 逐步筛选法（both）：综合了以上两种，直到AIC值最优。

什么是AIC值？
AIC又叫最小信息准则，是个日本人创立的，公式如下，总之，值越小越好，这里就不展开讲了，可以看下wiki里的介绍（点我）

$AIC = n \cdot ln(RSS_p / n) + 2p$
其中， $n$为变量总个数，$RSS$为残差平方和，p为选出的变量个数。
令人可喜的是，R语言中封装好了逐步回归的函数，我们直接调用即可！

s1 = step(lm.solMul2, direction = "both")

得到

我们来看下究竟发生了什么。刚开始的AIC值为151.74，如果剔除 $isEnglishYes*log(UniqueVisitors)$的话就变成了149.78，如果剔除 $hasAdYes*log(UniqueVisitors)$的话就变成了150.24，什么都不干还是151.74，于是自动剔除了 $isEnglishYes*log(UniqueVisitors)$。接着剔除 $hasAdYes*log(UniqueVisitors)$的话可以变成更小 148.28，其它方法只会使得AIC变大，故又剔除了 $hasAdYes*log(UniqueVisitors)$。最后，发现不管增加还是剔除变量都不能使得AIC减小了，于是逐步回归停止，我们得到了一个模型
$lny = a + b \cdot lnx + c \cdot IsEnglish + d \cdot HasAdvertising$
看下这个模型的效果怎么样~

summary(s1)

得到

怎么样~是不是挺不错的！每个参数的估计值显著性都很高，残差的标准差也挺小的，就是相关系数的平方值稍微低了那么一些，但还是可以令人接受的。好了，就要它了~
我们拿这个模型来看看预测的效果怎么样

predictData <- data.frame(UniqueVisitors = 280000000, hasAdYes = 1, isEnglishYes = 0)
lm.pred <- predict(s1, predictData, interval = "prediction", level = 0.95)
lm.pred

得到

也就是说我们认为结果是23.52291，并且结果落在(21.38347, 25.66234)这个区间的可能性为95%。
这里用的这组数据的真是值是

我们把PageViews的4.4e+10取个对数值得到

可见和我们的预测结果是差不多的~
当然，如果是在没取对数的情况下，这两者之间还是差的很远很远的~
所以，用多元线性回归预测网页访问量是不太准确的，可能考虑了其它的因素，寻找一下其它的特征之后再来算可能会好一点（可以试一下），换成其它的方法也可能会好些。
这里主要讲的的多元线性回归的解题思路，也就不管结果能不能用于实际情况了~

结束语

以上给出了解决多元线性回归问题的一般思路，该方法当然也适用于其他的案例。
如有不足，还请指正！

参考文献

[1] Drew Conway. 机器学习:实用案例解析[M]. 机械工业出版社, 2013.
[2] 薛毅. 统计建模与R软件[M]. 清华大学出版社, 2007.