R語言線性回歸預(yù)測(cè)網(wǎng)頁流量

 回歸是用已知的數(shù)據(jù)集來預(yù)測(cè)另一個(gè)數(shù)據(jù)集，如保險(xiǎn)精算師也許想在已知人們吸煙習(xí)慣的基礎(chǔ)上預(yù)測(cè)其壽命?；貧w模型的輸出是數(shù)字。
1、基準(zhǔn)模型
如果我們要在不使用其他任何信息的情況下，盡可能做出接近事實(shí)的預(yù)測(cè)，那么平均輸出作為結(jié)果是我們可以做的最好預(yù)測(cè)。在保險(xiǎn)精算師的例子中，我們可以完全忽略一個(gè)人的健康記錄并且預(yù)測(cè)其壽命等于人類平均壽命。

在討論如何做出最好的合理預(yù)測(cè)之前，假如我們有一組虛構(gòu)的保險(xiǎn)統(tǒng)計(jì)數(shù)據(jù)，第一列為是否抽煙（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

從這個(gè)圖中可以看出，吸煙習(xí)慣和壽命的關(guān)系，因?yàn)椴晃鼰煹膲勖植贾行暮臀鼰煹娜讼啾龋蛴移啤?br />
如果使用平方誤差作為預(yù)測(cè)質(zhì)量的衡量指標(biāo)，那么對(duì)人的壽命做出的最好假設(shè)（在沒有任何關(guān)于人的習(xí)慣信息的情況下）就是人的壽命均值。平方誤差的計(jì)算：（y-h）^2，其中y是真實(shí)結(jié)果，h是預(yù)測(cè)的結(jié)果。下面我們可以驗(yàn)證一下。

人的評(píng)價(jià)年齡可以使用mean方法獲得，在這里我們得到了72.723，向上取整得到73

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

通過假設(shè)已有數(shù)據(jù)集中每個(gè)人的壽命都是73，而得到的均方誤差是32.991。為了證明73是最好的假設(shè)，我們?cè)诜秶?3-83的可能假設(shè)序列上做一個(gè)循環(huán)。

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之外的其他任何假設(shè)對(duì)于我們的數(shù)據(jù)集來說帶來的都是更差的預(yù)測(cè)。這實(shí)際上是一個(gè)我們可以從數(shù)學(xué)上證明的一般理論結(jié)果：為了最小化均方誤差，需要使用數(shù)據(jù)集的均值作為預(yù)測(cè)。這說明了很重要的一點(diǎn)：在已有了關(guān)于吸煙信息的情況下做出預(yù)測(cè)，如果要衡量其好壞，那就應(yīng)該看它比你對(duì)每個(gè)人都用均值去猜的結(jié)果提升了多少。
2、使用虛擬變量的回歸模型

如何使用是否吸煙這樣的信息來對(duì)人的壽命做出更好的假設(shè)？一個(gè)簡單的想法是，先分別估算吸煙的人和不吸煙的人的死亡年齡均值，然后根據(jù)要研究的人是否吸煙，以對(duì)應(yīng)均值作為其預(yù)測(cè)壽命。這一次，我們使用均方根誤差（Root Mean Squared Error，RMSE）來代替均方誤差（MSE）。

下面是將吸煙的人和不吸煙的人分成單獨(dú)建模的兩組之后，使用R語言計(jì)算均方根誤差。

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  # 不包含吸煙信息的預(yù)測(cè)誤差

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  # 包含吸煙信息的預(yù)測(cè)誤差

從上例可以看出，在引入了更多的信息之后，所做出的預(yù)測(cè)確實(shí)更好了：當(dāng)引入關(guān)于吸煙習(xí)慣的信息之后，在預(yù)測(cè)人群壽命時(shí)的預(yù)測(cè)誤差減少了10%。

一般來說，每當(dāng)我們有了可以將數(shù)據(jù)點(diǎn)分為兩種類型的二元區(qū)分性質(zhì)–假設(shè)這些二元區(qū)分性和我們嘗試預(yù)測(cè)的結(jié)果相關(guān)，我們都能得到比僅僅使用均值更好的預(yù)測(cè)結(jié)果。簡單二元區(qū)分性的例子有：男人和女人。
3、線性回歸簡介

當(dāng)使用線性回歸模型預(yù)測(cè)輸出結(jié)果時(shí)，所做的最大的兩個(gè)假設(shè)如下：

1）可分性/可加性

如果有多份信息可能影響我們的假設(shè)，那么通過累加每一份信息的影響來產(chǎn)生我們的假設(shè)，就像單獨(dú)使用每份信息時(shí)一樣。假如，如果酗酒者比不酗酒者少活1年，并且吸煙者比不吸煙者少活5年，那么一個(gè)吸煙的酗酒者應(yīng)該會(huì)比既不吸煙也不酗酒的人少活6（1+5）年。這種假設(shè)是事情同時(shí)發(fā)生時(shí)將他們的單獨(dú)影響累加在一起，是一個(gè)很大的假設(shè)，但是這是很多回歸模型應(yīng)用的不錯(cuò)起點(diǎn)。

2）單調(diào)性/線性

當(dāng)改變一個(gè)輸入值總使得預(yù)測(cè)的結(jié)果增加或者減少時(shí)，這個(gè)模型是單調(diào)的。假如，你使用身高作為輸入值預(yù)測(cè)體重，并且模型是單調(diào)的，那么當(dāng)前的預(yù)測(cè)是每當(dāng)某些人的身高增加，他們的體重將會(huì)增加。如果將輸入和輸出畫出來，將會(huì)看到一條直線，而不是某種更復(fù)雜的形狀，如曲線或者波浪線。

使用身高體重的例子，在調(diào)用geom_smooth函數(shù)時(shí)指明要使用lm方法即可，其中l(wèi)m方法已經(jīng)實(shí)現(xiàn)了“線性模型”。

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

從上圖中可以看出，通過這條直線，在已知一個(gè)人身高的前提下去預(yù)測(cè)其體重回去的非常好的效果。例如，看著這條直線，我們可以預(yù)測(cè)身高60英寸的人體重為105磅。至于如何找到用于定義在這幅畫種看到的直線的數(shù)字，這正是R語言所擅長的地方：R語言中一個(gè)稱為lm的簡單函數(shù)將會(huì)為我們完成所有的這些工作。為了使用lm，我們需要使用~操作符指明一個(gè)公式。

我們可以使用下面的公式運(yùn)行一個(gè)線性回歸程序：

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

一旦運(yùn)行了lm函數(shù)的調(diào)用，就可以通過調(diào)用coef函數(shù)來得到回歸直線的截距。

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

這也就是說某個(gè)人的身高增加一英寸，就會(huì)導(dǎo)致他的體重增加7.7磅。但是這個(gè)模型中，一個(gè)人至少有45英寸身高，才能顯示出其體重0磅。簡言之，我們的回歸模型對(duì)于兒童或者身高特別矮的成年人來說并不是太使用。

predict(fitted.regression)

predict可以獲得模型對(duì)于每個(gè)數(shù)值的預(yù)測(cè)結(jié)果。一旦有了這個(gè)結(jié)果，就可以使用簡單的減法來計(jì)算預(yù)測(cè)結(jié)果和真實(shí)值之間的誤差。

true.values <- with(heights.weights, Weight)
errors <- true.values - predict(fitted.regression)  # 真實(shí)值和預(yù)測(cè)值之間的差，也叫殘差

在R語言中可以使用residuals函數(shù)替代predict函數(shù)來直接獲得殘差：

> 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

為了發(fā)現(xiàn)使用線性回歸時(shí)產(chǎn)生的明顯錯(cuò)誤，可以把殘差和真實(shí)數(shù)據(jù)對(duì)應(yīng)畫在一幅畫中。

plot(fitted.regression, which = 1)
#which=1盡讓R語言畫出了第一個(gè)回歸診斷點(diǎn)圖。

3

在這個(gè)例子中，我們可以說這個(gè)線性模型很有效，因?yàn)闅埐钪胁淮嬖谙到y(tǒng)性的結(jié)構(gòu)（？？沒看明白）。

下面舉一個(gè)直線并不適用的例子：

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

4

對(duì)于這個(gè)問題，我們可以看到殘差中存在明顯的結(jié)構(gòu)。

最簡單的誤差衡量指標(biāo)是：1）取得所有的殘差；2）對(duì)他們進(jìn)行平方處理，以獲取模型的誤差平方；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

對(duì)于比較不同的模型，這個(gè)簡單的誤差平方和數(shù)值是有用的，但是誤差平方和在大數(shù)據(jù)集上的值比在小數(shù)據(jù)集上的值更大。我們可以使用誤差平方的均值來代替這個(gè)誤差平方和，也就是前面提到過的均方誤差（MSE）度量方法。

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

對(duì)均方誤差進(jìn)行開放運(yùn)算以獲得均方根誤差，這就是RMES度量方法，這宗方法一般用于評(píng)估機(jī)器學(xué)習(xí)算法的效果。
RMSE有一點(diǎn)不盡人意，就是它不能讓人直觀清晰地看出哪個(gè)模型表現(xiàn)平平。理想的效果是RMSE值為0。同樣，使用RMSE也不容易識(shí)別什么時(shí)候一個(gè)模型的效果非常差。例如，如果每個(gè)人的身高都是5英寸，而預(yù)測(cè)結(jié)果是5000英寸，這時(shí)見得到一個(gè)巨大的RMSE。為了解決這個(gè)問題可以使用R2，下面說明了如何得到R2，第一步只是用均值來當(dāng)做所有樣本數(shù)據(jù)的預(yù)測(cè)值時(shí)的RMSE，第二步是使用你的模型所作出的預(yù)測(cè)的RMSE。

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

4、預(yù)測(cè)網(wǎng)頁流量
在這里我們使用線性回歸模型預(yù)測(cè)互聯(lián)網(wǎng)上排名前1000的網(wǎng)站在2011年的訪問量。數(shù)據(jù)項(xiàng)主要分布如下：

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是網(wǎng)站的排名，PageViewss是一年中網(wǎng)站被訪問了多少次，UniqueVisitors是有多少不同的用戶訪問網(wǎng)站，HasAdVertising一個(gè)網(wǎng)站上是否有廣告，IsEnglish是網(wǎng)站上的語言是否為英語。

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

5
從上圖中可以看到，幾乎所有的數(shù)據(jù)都在X軸的附近集成一束，這是使用非標(biāo)準(zhǔn)分布數(shù)據(jù)工作時(shí)常見的一個(gè)問題。下面我們來看看PageViews本身的分布：

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

6
這個(gè)密度圖和前面的散點(diǎn)圖一樣不可理解，當(dāng)看到?jīng)]有意義的密度圖時(shí)，最好的方法是嘗試對(duì)你想要分析的數(shù)值取log，并且經(jīng)過log后重新繪制一幅密度圖。
7
這樣的密度圖看起來就合理多了，因此我們就使用log變換后的PageView和UniqueVisitors。散點(diǎn)圖的作圖結(jié)果如下圖所示，看上去好像有一條可以使用回歸模型畫出的潛在的直線。我們以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函數(shù)來找到定義這條直線斜率和截距的數(shù)值：

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函數(shù)告訴我們的第一件事情是對(duì)lm函數(shù)所做的調(diào)用。當(dāng)使用對(duì)lm進(jìn)行了多次調(diào)用的大型腳本進(jìn)行工作時(shí)，該參數(shù)就變得非常有用。

summary函數(shù)告訴我們的第二件事情是殘差的分位數(shù)。如果調(diào)用quantile(residuals(lm.fit))也可以計(jì)算出這個(gè)分位數(shù)。

接著，summary提供了比coef函數(shù)更詳細(xì)的回歸模型系數(shù)信息。每一個(gè)系數(shù)都有一個(gè)Estimate，Std. Error， t-value， Pr(>|t|)。這些值用于評(píng)估我們計(jì)算結(jié)果存在的不確定性，換句話說，他們是置信度。如“Std.Error”可以用于產(chǎn)生一個(gè)置信度為95%的系數(shù)置信區(qū)間?！皌-value”和“p-value”用于衡量我們對(duì)真實(shí)系數(shù)不為零有多大信心。在本例中，log(UniqueVisitors)的系數(shù)是1.33628，而標(biāo)準(zhǔn)差是0.04568，就是說這個(gè)系數(shù)距離零26.25306（1.33628/0.04568 == 26.25306）。如果得到的系數(shù)與零距離遠(yuǎn)在3個(gè)標(biāo)簽誤差之上，那么就有理由相信這兩個(gè)變量之間是相關(guān)的。

下一部分信息是關(guān)于系數(shù)的顯著性編碼。數(shù)字旁邊的星號(hào)的意思是“t-value”有多大或者“p-value”有多小。

最后一部分信息室關(guān)于從數(shù)據(jù)中擬合得到的現(xiàn)行模型的預(yù)測(cè)能力。第一個(gè)是“Residual standard error”，就是使用sqrt(mean(residuals(lm.fit)^2))計(jì)算出來的RMSE?！癲egrees of freedom”我們?cè)诜治鲋惺褂玫臄?shù)據(jù)點(diǎn)至少要有兩個(gè)，才能有效地?cái)M合兩個(gè)系數(shù)。“Multiple R-squared”是標(biāo)準(zhǔn)的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