R語言線性回歸預測網頁流量

 回歸是用已知的數據集來預測另一個數據集，如保險精算師也許想在已知人們吸煙習慣的基礎上預測其壽命?；貧w模型的輸出是數字。
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%的系數置信區間?！皌-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?！癲egrees of freedom”我們在分析中使用的數據點至少要有兩個，才能有效地擬合兩個系數?！癕ultiple 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