R语言线性回归预测网页流量

 回归是用已知的数据集来预测另一个数据集，如保险精算师也许想在已知人们吸烟习惯的基础上预测其寿命。回归模型的输出是数字。
1、基准模型
如果我们要在不使用其他任何信息的情况下，尽可能做出接近事实的预测，那么平均输出作为结果是我们可以做的最好预测。在保险精算师的例子中，我们可以完全忽略一个人的健康记录并且预测其寿命等于人类平均寿命。

在讨论如何做出最好的合理预测之前，假如我们有一组虚构的保险统计数据，第一列为是否抽烟（0不抽烟，1为抽烟），第二列是年龄。我们先用密度图来比较吸烟者和非吸烟者，如下所示。

library('ggplot2')
# First snippet
ages <- read.csv(file.path('data', 'longevity.csv'))
ggplot(ages, aes(x = AgeAtDeath, fill = factor(Smokes))) +
  geom_density() +
  facet_grid(Smokes ~ .)

g

从这个图中可以看出，吸烟习惯和寿命的关系，因为不吸烟的寿命分布中心和吸烟的人相比，向右偏移。

如果使用平方误差作为预测质量的衡量指标，那么对人的寿命做出的最好假设（在没有任何关于人的习惯信息的情况下）就是人的寿命均值。平方误差的计算：（y-h）^2，其中y是真实结果，h是预测的结果。下面我们可以验证一下。

人的评价年龄可以使用mean方法获得，在这里我们得到了72.723，向上取整得到73

> mean(ages$AgeAtDeath)
[1] 72.723
> guess <- 73
> with(ages,mean((AgeAtDeath - guess) ^2))
[1] 32.991

通过假设已有数据集中每个人的寿命都是73，而得到的均方误差是32.991。为了证明73是最好的假设，我们在范围63-83的可能假设序列上做一个循环。

guess.accuracy <- data.frame()

for (guess in seq(63, 83, by = 1))
{
 prediction.error <- with(ages,
 mean((AgeAtDeath - guess) ^ 2))
 guess.accuracy <- rbind(guess.accuracy,
 data.frame(Guess = guess,
 Error = prediction.error))
}

ggplot(guess.accuracy, aes(x = Guess, y = Error)) +
 geom_point() +
 geom_line()

1

如上图所示，使用除了73之外的其他任何假设对于我们的数据集来说带来的都是更差的预测。这实际上是一个我们可以从数学上证明的一般理论结果：为了最小化均方误差，需要使用数据集的均值作为预测。这说明了很重要的一点：在已有了关于吸烟信息的情况下做出预测，如果要衡量其好坏，那就应该看它比你对每个人都用均值去猜的结果提升了多少。
2、使用虚拟变量的回归模型

如何使用是否吸烟这样的信息来对人的寿命做出更好的假设？一个简单的想法是，先分别估算吸烟的人和不吸烟的人的死亡年龄均值，然后根据要研究的人是否吸烟，以对应均值作为其预测寿命。这一次，我们使用均方根误差（Root Mean Squared Error，RMSE）来代替均方误差（MSE）。

下面是将吸烟的人和不吸烟的人分成单独建模的两组之后，使用R语言计算均方根误差。

ages <- read.csv(file.path('data', 'longevity.csv'))
constant.guess <- with(ages, mean(AgeAtDeath))
with(ages, sqrt(mean((AgeAtDeath - constant.guess) ^ 2)))
# [1] 5.737096  # 不包含吸烟信息的预测误差

smokers.guess <- with(subset(ages, Smokes == 1),mean(AgeAtDeath))
non.smokers.guess <- with(subset(ages, Smokes == 0),
                          mean(AgeAtDeath))
ages <- transform(ages,
                  NewPrediction = ifelse(Smokes == 0,
                  non.smokers.guess,
                  smokers.guess))
with(ages, sqrt(mean((AgeAtDeath - NewPrediction) ^ 2)))
# [1] 5.148622  # 包含吸烟信息的预测误差

从上例可以看出，在引入了更多的信息之后，所做出的预测确实更好了：当引入关于吸烟习惯的信息之后，在预测人群寿命时的预测误差减少了10%。

一般来说，每当我们有了可以将数据点分为两种类型的二元区分性质–假设这些二元区分性和我们尝试预测的结果相关，我们都能得到比仅仅使用均值更好的预测结果。简单二元区分性的例子有：男人和女人。
3、线性回归简介

当使用线性回归模型预测输出结果时，所做的最大的两个假设如下：

1）可分性/可加性

如果有多份信息可能影响我们的假设，那么通过累加每一份信息的影响来产生我们的假设，就像单独使用每份信息时一样。假如，如果酗酒者比不酗酒者少活1年，并且吸烟者比不吸烟者少活5年，那么一个吸烟的酗酒者应该会比既不吸烟也不酗酒的人少活6（1+5）年。这种假设是事情同时发生时将他们的单独影响累加在一起，是一个很大的假设，但是这是很多回归模型应用的不错起点。

2）单调性/线性

当改变一个输入值总使得预测的结果增加或者减少时，这个模型是单调的。假如，你使用身高作为输入值预测体重，并且模型是单调的，那么当前的预测是每当某些人的身高增加，他们的体重将会增加。如果将输入和输出画出来，将会看到一条直线，而不是某种更复杂的形状，如曲线或者波浪线。

使用身高体重的例子，在调用geom_smooth函数时指明要使用lm方法即可，其中lm方法已经实现了“线性模型”。

library('ggplot2')
heights.weights <- read.csv(file.path('data',
                                      '01_heights_weights_genders.csv'),
                             header = TRUE,
                             sep = ',')
ggplot(heights.weights, aes(x = Height, y = Weight)) +
  geom_point() +
  geom_smooth(method = 'lm')

2

从上图中可以看出，通过这条直线，在已知一个人身高的前提下去预测其体重回去的非常好的效果。例如，看着这条直线，我们可以预测身高60英寸的人体重为105磅。至于如何找到用于定义在这幅画种看到的直线的数字，这正是R语言所擅长的地方：R语言中一个称为lm的简单函数将会为我们完成所有的这些工作。为了使用lm，我们需要使用~操作符指明一个公式。

我们可以使用下面的公式运行一个线性回归程序：

fitted.regression <- lm(Weight ~ Height,data = heights.weights)

一旦运行了lm函数的调用，就可以通过调用coef函数来得到回归直线的截距。

coef(fitted.regression)
#(Intercept) Height
#-350.737192 7.717288
# predicted.weight == -350.737192 + 7.717288 * observed.height

这也就是说某个人的身高增加一英寸，就会导致他的体重增加7.7磅。但是这个模型中，一个人至少有45英寸身高，才能显示出其体重0磅。简言之，我们的回归模型对于儿童或者身高特别矮的成年人来说并不是太使用。

predict(fitted.regression)

predict可以获得模型对于每个数值的预测结果。一旦有了这个结果，就可以使用简单的减法来计算预测结果和真实值之间的误差。

true.values <- with(heights.weights, Weight)
errors <- true.values - predict(fitted.regression)  # 真实值和预测值之间的差，也叫残差

在R语言中可以使用residuals函数替代predict函数来直接获得残差：

> head(heights.weights)
 Gender Height Weight
1 Male 73.84702 241.8936
2 Male 68.78190 162.3105
3 Male 74.11011 212.7409
4 Male 71.73098 220.0425
5 Male 69.88180 206.3498
6 Male 67.25302 152.2122
> head(predict(fitted.regression))
 1 2 3 4 5 6
219.1615 180.0725 221.1918 202.8314 188.5607 168.2737
> head(errors)
 1 2 3 4 5 6
 22.732083 -17.762074 -8.450953 17.211069 17.789073 -16.061519
> head(residuals(fitted.regression))
 1 2 3 4 5 6
 22.732083 -17.762074 -8.450953 17.211069 17.789073 -16.061519

为了发现使用线性回归时产生的明显错误，可以把残差和真实数据对应画在一幅画中。

plot(fitted.regression, which = 1)
#which=1尽让R语言画出了第一个回归诊断点图。

3

在这个例子中，我们可以说这个线性模型很有效，因为残差中不存在系统性的结构（？？没看明白）。

下面举一个直线并不适用的例子：

x <- 1:10
y <- x ^ 2
fitted.regression <- lm(y ~ x)
plot(fitted.regression, which = 1)

4

对于这个问题，我们可以看到残差中存在明显的结构。

最简单的误差衡量指标是：1）取得所有的残差；2）对他们进行平方处理，以获取模型的误差平方；3）把这些误差平方加载一起求和。

x <- 1:10
y <- x ^ 2
fitted.regression <- lm(y ~ x)
errors <- residuals(fitted.regression)
squared.errors <- errors ^ 2
sum(squared.errors)
#[1] 528

对于比较不同的模型，这个简单的误差平方和数值是有用的，但是误差平方和在大数据集上的值比在小数据集上的值更大。我们可以使用误差平方的均值来代替这个误差平方和，也就是前面提到过的均方误差（MSE）度量方法。

mse <- mean(squared.errors)
mse
#[1] 52.8
rmse <- sqrt(mse)
rmse
#[1] 7.266361

对均方误差进行开放运算以获得均方根误差，这就是RMES度量方法，这宗方法一般用于评估机器学习算法的效果。
RMSE有一点不尽人意，就是它不能让人直观清晰地看出哪个模型表现平平。理想的效果是RMSE值为0。同样，使用RMSE也不容易识别什么时候一个模型的效果非常差。例如，如果每个人的身高都是5英寸，而预测结果是5000英寸，这时见得到一个巨大的RMSE。为了解决这个问题可以使用R2，下面说明了如何得到R2，第一步只是用均值来当做所有样本数据的预测值时的RMSE，第二步是使用你的模型所作出的预测的RMSE。

mean.mse <- 1.09209343
model.mse <- 0.954544
r2 <- 1 - (model.mse / mean.mse)
r2
#[1] 0.1259502

4、预测网页流量
在这里我们使用线性回归模型预测互联网上排名前1000的网站在2011年的访问量。数据项主要分布如下：

top.1000.sites <- read.csv(file.path('data', 'top_1000_sites.tsv'),
                           sep = '\t',
                           stringsAsFactors = FALSE)

 >head(top.1000.sites)
 Rank Site Category UniqueVisitors Reach PageViews HasAdvertising InEnglish TLD
1 1 facebook.com   Social Networks  880000000 47.2 9.1e+11 Yes Yes com
2 2 youtube.com    Online Video     800000000 42.7 1.0e+11 Yes Yes com
3 3 yahoo.com      Web Portals      660000000 35.3 7.7e+10 Yes Yes com
4 4 live.com       Search Engines   550000000 29.3 3.6e+10 Yes Yes com
5 5 wikipedia.org  Dictionaries & Encyclopedias  490000000 26.2 7.0e+09 No Yes org
6 6 msn.com        Web Portals      450000000 24.0 1.5e+10 Yes Yes com

我们主要考虑如下的五列：Rank、PageViews、UniqueVisitors、HasAdVertising和IsEnglish。其中Rank是网站的排名，PageViewss是一年中网站被访问了多少次，UniqueVisitors是有多少不同的用户访问网站，HasAdVertising一个网站上是否有广告，IsEnglish是网站上的语言是否为英语。

ggplot(top.1000.sites, aes(x = PageViews, y = UniqueVisitors)) +
  geom_point()

5
从上图中可以看到，几乎所有的数据都在X轴的附近集成一束，这是使用非标准分布数据工作时常见的一个问题。下面我们来看看PageViews本身的分布：

ggplot(top.1000.sites, aes(x = PageViews)) +geom_density()

6
这个密度图和前面的散点图一样不可理解，当看到没有意义的密度图时，最好的方法是尝试对你想要分析的数值取log，并且经过log后重新绘制一幅密度图。
7
这样的密度图看起来就合理多了，因此我们就使用log变换后的PageView和UniqueVisitors。散点图的作图结果如下图所示，看上去好像有一条可以使用回归模型画出的潜在的直线。我们以method=‘lm’使用geom_smooth来看看回归直线将是什么样子的：

ggplot(top.1000.sites, aes(x = log(PageViews), y = log(UniqueVisitors))) +
  geom_point() +
  geom_smooth(method = 'lm', se = FALSE)

8

我们可以通过lm函数来找到定义这条直线斜率和截距的数值：

lm.fit <- lm(log(PageViews) ~ log(UniqueVisitors),data = top.1000.sites)
# Twenty-third snippet
summary(lm.fit)
#Call:
#lm(formula = log(PageViews) ~ log(UniqueVisitors), data = top.1000.sites)
#
#Residuals:
# Min 1Q Median 3Q Max
#-2.1825 -0.7986 -0.0741 0.6467 5.1549
#
#Coefficients:
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) -2.83441 0.75201 -3.769 0.000173 ***
#log(UniqueVisitors) 1.33628 0.04568 29.251 < 2e-16 ***
#---
#Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#Residual standard error: 1.084 on 998 degrees of freedom
#Multiple R-squared: 0.4616, Adjusted R-squared: 0.4611
#F-statistic: 855.6 on 1 and 998 DF, p-value: < 2.2e-16

summary函数告诉我们的第一件事情是对lm函数所做的调用。当使用对lm进行了多次调用的大型脚本进行工作时，该参数就变得非常有用。

summary函数告诉我们的第二件事情是残差的分位数。如果调用quantile(residuals(lm.fit))也可以计算出这个分位数。

接着，summary提供了比coef函数更详细的回归模型系数信息。每一个系数都有一个Estimate，Std. Error， t-value， Pr(>|t|)。这些值用于评估我们计算结果存在的不确定性，换句话说，他们是置信度。如“Std.Error”可以用于产生一个置信度为95%的系数置信区间。“t-value”和“p-value”用于衡量我们对真实系数不为零有多大信心。在本例中，log(UniqueVisitors)的系数是1.33628，而标准差是0.04568，就是说这个系数距离零26.25306（1.33628/0.04568 == 26.25306）。如果得到的系数与零距离远在3个标签误差之上，那么就有理由相信这两个变量之间是相关的。

下一部分信息是关于系数的显著性编码。数字旁边的星号的意思是“t-value”有多大或者“p-value”有多小。

最后一部分信息室关于从数据中拟合得到的现行模型的预测能力。第一个是“Residual standard error”，就是使用sqrt(mean(residuals(lm.fit)^2))计算出来的RMSE。“degrees of freedom”我们在分析中使用的数据点至少要有两个，才能有效地拟合两个系数。“Multiple R-squared”是标准的R平方。

# Twenty-fourth snippet
lm.fit <- lm(log(PageViews) ~ HasAdvertising + log(UniqueVisitors) + InEnglish,
 data = top.1000.sites)

summary(lm.fit)
#Call:
#lm(formula = log(PageViews) ~ HasAdvertising + log(UniqueVisitors) +
# InEnglish, data = top.1000.sites)
#
#Residuals:
# Min 1Q Median 3Q Max
#-2.4283 -0.7685 -0.0632 0.6298 5.4133
#
#Coefficients:
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) -1.94502 1.14777 -1.695 0.09046 .
#HasAdvertisingYes 0.30595 0.09170 3.336 0.00088 ***
#log(UniqueVisitors) 1.26507 0.07053 17.936 < 2e-16 ***
#InEnglishNo 0.83468 0.20860 4.001 6.77e-05 ***
#InEnglishYes -0.16913 0.20424 -0.828 0.40780
#---
#Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#Residual standard error: 1.067 on 995 degrees of freedom
#Multiple R-squared: 0.4798, Adjusted R-squared: 0.4777
#F-statistic: 229.4 on 4 and 995 DF, p-value: < 2.2e-16