1.線性回歸入門

線性回歸的假設(shè)

假設(shè)輸出變量是一組特征變量的加權(quán)線性函數(shù)

?假設(shè)1：誤差項(xiàng)具有同方差性

同方差（homoscedasticity）：誤差部分的方差不會(huì)隨著輸入特征的值或水平而變化
異方差（heteroskedastic）：誤差部分的方差隨著輸入特征的值的變化而變化

order()函數(shù)用于返回向量大小順序的秩溜腐，從小到大的數(shù)字的序號(hào)值调俘。
x<-c(10,6,4,7,8)
order(x)
[1] 3 2 4 5 1
sapply(x,FUN)臀稚，第一個(gè)參數(shù)x是需要處理的數(shù)據(jù)疫粥，F(xiàn)UN是處理數(shù)據(jù)x的Funtion蟆湖。

library(ggplot2)
set.seed(15)
# 隨機(jī)產(chǎn)生200個(gè)數(shù)故爵，均值100，標(biāo)準(zhǔn)差25
x1 <- order(rnorm(200, mean = 100, sd = 25))
y <- sapply(x1, function(x) {17+1.8*x + rnorm(1,mean=0, sd = x/4)})
#標(biāo)準(zhǔn)差隨x變大而變大隅津，y的值也隨之變大诬垂，誤差項(xiàng)的方差與x有相關(guān)性劲室。其實(shí)是構(gòu)造了一個(gè)隨x變化的異方差的y的數(shù)據(jù)
df <- data.frame(x1 = x1, y = y)
fit<-coef(lm(y ~ x1, data = df))
# coef()函數(shù)提取線性回歸lm模型的系數(shù)
# qplot畫點(diǎn)圖，geom_abline畫線圖
# size：坐標(biāo)軸標(biāo)簽字體大小结窘，lineheight: 標(biāo)簽行間距的倍數(shù)很洋，family:字體，face:字體外形（粗斜體等）
p <- qplot(x1,y,data=df) + geom_abline(intercept=fit[1],slope=fit[2]) + ggtitle("Linear Relationship with Heteroskedastic Errors")+      
theme(plot.title = element_text(lineheight=.8, face="bold"))
p

可以看到隧枫，系數(shù)是17.0848喉磁，1.81745，非常接近我們構(gòu)造的函數(shù)17+1.8*x + rnorm(1,mean=0, sd = x/4)悠垛，因?yàn)橛性肼曄叨ǎ越鼐嗪蛒的系數(shù)有細(xì)微變化。

同方差性與異方差性的區(qū)別：

a <- 1:200
#同方差
A <- data.frame(x=1:200, y=1+2*a+rnorm(n=200,mean=0,sd=10))
ggplot(data = A, aes(x,y)) +
  geom_point()+ geom_smooth(method = 'lm')
#異方差
B <- data.frame(x=1:200, y=1+2*a+rnorm(n=200,mean=0,sd=0.4*a))
ggplot(data = B, aes(x,y)) +
  geom_point()+ geom_smooth(method = 'lm')

?假設(shè)2：誤差項(xiàng)具有獨(dú)立性

人工擴(kuò)大數(shù)據(jù)集的方式會(huì)破壞獨(dú)立性的假設(shè)确买。
例如：取出一部分?jǐn)?shù)據(jù)斤讥，把特征值和輸出值都乘以2，看起來(lái)數(shù)據(jù)集合變大了湾趾，但是數(shù)據(jù)之間的誤差有了相關(guān)性芭商。
不同工廠的測(cè)量數(shù)據(jù)，工廠應(yīng)該作為一種特征搀缠。

?假設(shè)3：特征在統(tǒng)計(jì)學(xué)上是相互獨(dú)立的

2.簡(jiǎn)單線性回歸

從均勻分布中隨機(jī)抽取點(diǎn)铛楣，讓數(shù)據(jù)理想地分開
?? = 1.67??1 ? 2.93 + ??(0,2^2)

估算回歸系數(shù)

注意：
如果 ??1^ = 0表示輸出變量與輸入特征相互獨(dú)立
即使兩個(gè)變量統(tǒng)計(jì)學(xué)上相互獨(dú)立，通常還是會(huì)有少量的協(xié)方差
協(xié)方差是0也不代表是相互獨(dú)立
如果訓(xùn)練一個(gè)線性回歸模型艺普，一般??1^ ≠ 0

set.seed(5427395)
nObs <- 100
x1minrange <- 5 
x1maxrange <- 25
# 隨機(jī)產(chǎn)生100個(gè)最小值為5簸州，最大值為25的均勻分布隨機(jī)數(shù)
x1 <- runif(nObs,x1minrange,x1maxrange)
# 隨機(jī)產(chǎn)生100個(gè)均值為0，標(biāo)準(zhǔn)差為2的正態(tài)分布隨機(jī)數(shù)
e <- rnorm(nObs,mean = 0, sd = 2.0)
y <- 1.67*x1 - 2.93 + e
# 構(gòu)造一個(gè)數(shù)據(jù)框格式
df <- data.frame(y,x1)
# 線性擬合
myfit <- lm(y~x1,df)
summary(myfit)

因?yàn)閑的噪聲干擾歧譬，發(fā)現(xiàn)截距在-2.93附近岸浑，x1自變量系數(shù)在1.67附近。e是同方差干擾瑰步。

# 畫圖
p <- qplot(x1,y,data=df) 
p <- p + ggtitle("Simple Linear Regression")  #畫標(biāo)題
# 畫lm函數(shù)擬合直線
p <- p + geom_abline(intercept=coef(myfit)[1],slope=coef(myfit)[2], aes(linetype="Estimated Line"), size=1.2, show_guide=T)  
# 畫y = 1.67*x1 - 2.93 + e固定系數(shù)下原始直線
p <- p + geom_abline(intercept = -2.93, slope = 1.67, aes(linetype="Population Regression Line"),size=1.2, show_guide=T) 
p <- p + xlab("Feature x1")  # x軸的名稱
p <- p + ylab("Output y")  # y軸的名稱
# 控制線性矢洲，也可以寫scale_linetype_manual(values=c("twodash", "dotted"))
p <- p + scale_linetype_manual(name="",values=c("Population Regression Line"=1,"Estimated Line"=2))  
p <- p + theme(plot.title = element_text(lineheight=.8, face="bold"),legend.position = "bottom")
p

以上是回歸線性系數(shù)估計(jì)擬合線與固定系數(shù)和截距劃線的比較，發(fā)現(xiàn)兩條線基本重合缩焦，擬合度較好读虏。

3.多元線性回歸

具有k個(gè)輸入特征的多元線性回歸模型?? = ???????? + ?????1?????1 + ? + ??1??1 + ??0 + ??
通常假設(shè)??與(??1, ??2, … , ????)相互獨(dú)立。

1）預(yù)測(cè)CPU性能

把數(shù)據(jù)集做成訓(xùn)練集和測(cè)試集袁滥，比例85%:15%
1.該數(shù)據(jù)沒(méi)有缺失值盖桥，而且樣本較少，僅有200行

library(caret)
set.seed(4352345)
machine$PRP
# 按照輸出變量PRP進(jìn)行索引题翻，提取85%作為訓(xùn)練集
machine_sampling_vector <- createDataPartition(machine$PRP, p = 0.85, list = FALSE)
# 獲取訓(xùn)練集數(shù)據(jù)
machine_train <- machine[machine_sampling_vector,]
machine_train_labels <- machine$PRP[machine_sampling_vector]

# 獲取測(cè)試集數(shù)據(jù)
machine_test <- machine[-machine_sampling_vector,]
machine_test_labels <- machine$PRP[-machine_sampling_vector]

2.數(shù)據(jù)設(shè)置好后葱轩，分析和檢查針對(duì)線性回歸的假設(shè)

machine_train_features <-machine[,1:6] #保留1-6列，去掉了第七列PRP
machine_correlations <- cor(machine_train_features) #每一列之間的相關(guān)性

# 篩選與其他相關(guān)性系數(shù)強(qiáng)的變量并且刪除
findCorrelation(machine_correlations) 
findCorrelation(machine_correlations, cutoff = 0.75)
cor(machine_train$MMIN,machine_train$MMAX)

高度相關(guān)性的閾值默認(rèn)的0.9沒(méi)有發(fā)現(xiàn)哪個(gè)特征該去掉。
但是靴拱，當(dāng)閾值改成0.75的時(shí)候垃喊，caret包推薦去掉第三個(gè)特征MMAX。
檢查發(fā)現(xiàn)MMIN和MMAX相關(guān)性很大袜炕。最小主內(nèi)存的型號(hào)也有較大的主內(nèi)存本谜，這里我們不修改。

2）預(yù)測(cè)二手汽車的價(jià)格

數(shù)據(jù)集caret偎窘，獲取二手車可靠?jī)r(jià)格的數(shù)據(jù)乌助，包含804種通用汽車（GM）品牌的汽車。

第一步：數(shù)據(jù)預(yù)處理

library(caret)
data(cars)
# 查看一下各特征之間的相關(guān)系數(shù)矩陣陌知，判斷相關(guān)關(guān)系
cars_cor <- cor(cars[,-1]) #去掉第一列price
findCorrelation(cars_cor) #默認(rèn)cutoff=0.9
findCorrelation(cars_cor, cutoff=0.75)
cor(cars$Doors,cars$coupe)
# 用交叉列聯(lián)表來(lái)查看相關(guān)性
table(cars$coupe,cars$Doors)
# 查找完全線性組合的特征
findLinearCombos(cars)
# 根據(jù)建議他托，去掉具有完全線性組合的特征
cars <- cars[,c(-15,-18)] 

第三列doors有超過(guò)75%相關(guān)性。查看相關(guān)矩陣仆葡，發(fā)現(xiàn)在Doors和coupe之間存在較高的相關(guān)性赏参。
跑車大概率是2門車，跑車和車門相關(guān)性為-0.8254435沿盅。門越少把篓，是跑車可能性越高。2門車中140個(gè)是跑車腰涧，50個(gè)不是跑車韧掩，4門車中0個(gè)是跑車。

函數(shù)結(jié)果表示：要把15(coupe)和18(wagon)remove掉窖铡，因?yàn)?5和18是其他特征的組合疗锐。

第二步：把數(shù)據(jù)分為訓(xùn)練集和測(cè)試集

set.seed(232455)
# 劃分訓(xùn)練集和測(cè)試集焰宣，同時(shí)標(biāo)注特征數(shù)據(jù)和標(biāo)簽數(shù)據(jù)
cars_sampling_vector <- createDataPartition(cars$Price, p = 0.85, list = FALSE)
# 獲取訓(xùn)練集數(shù)據(jù)
cars_train <- cars[cars_sampling_vector,]
cars_train_features <- cars[,-1] #去掉第一列price
cars_train_labels <- cars$Price[cars_sampling_vector]
# 獲取測(cè)試集數(shù)據(jù)
cars_test <- cars[-cars_sampling_vector,]
cars_test_labels <- cars$Price[-cars_sampling_vector] 

4.評(píng)估線性回歸模型

利用lm函數(shù)畴椰，用線性回歸模型來(lái)擬合數(shù)據(jù)。

# 根據(jù)訓(xùn)練集建立模型
machine_model1 <- lm(PRP~.,data=machine_train)
cars_model1 <- lm(Price~.,data=cars_train)
# 模型評(píng)估（包括：殘差分析度迂、顯著性檢驗(yàn)等敌买。其中顯著性檢驗(yàn)又包括線性關(guān)系檢驗(yàn)和回歸系數(shù)檢驗(yàn)）
summary(cars_model1)
summary(machine_model1)

summary()結(jié)果解讀

1. 調(diào)用：Call，表明lm是如何被調(diào)用的
2. 殘差統(tǒng)計(jì)量：Residuals阶界，殘差第一四分位數(shù)（1Q）和第三分位數(shù)（Q3）有大約相同的幅度虹钮，意味著有較對(duì)稱的鐘形分布
3. 系數(shù)：Coefficients
   分別表示： 估值標(biāo)準(zhǔn)誤差 T值 P值
   Intercept：表示截距
   Mileage：影響因子/特征
   Estimate：包含由普通最小二乘法計(jì)算出來(lái)的估計(jì)回歸系數(shù)。
   Std. Error：估計(jì)的回歸系數(shù)的標(biāo)準(zhǔn)誤差膘融。
   t值：系數(shù)估計(jì)值與標(biāo)準(zhǔn)誤差的比值芙粱。
   P值：估計(jì)系數(shù)不顯著的可能性，有較大P值的變量是可以從模型中移除的候選變量氧映〈号希可以直接通過(guò)P值與預(yù)設(shè)的0.05進(jìn)行比較，來(lái)判定對(duì)應(yīng)的解釋變量的顯著性，檢驗(yàn)的原假設(shè)是：該系數(shù)顯著為0律姨；若P<0.05振峻，則拒絕原假設(shè)，即對(duì)應(yīng)的變量顯著不為0择份。（p值很高并不一定代表特征值和輸出沒(méi)有線性關(guān)系扣孟；它只表明有其他特征存在的時(shí)候，這個(gè)特征對(duì)輸出變量不提供新的信息）
   可以看到Mileage荣赶、Cylinder凤价、Doors都可以認(rèn)為是在P為0.05的水平下顯著不為0，通過(guò)顯著性檢驗(yàn)拔创；Cruise的P值為0.20025利诺，不顯著。
4. Multiple R-squared和Adjusted R-squared
   這兩個(gè)值剩燥，即R^{2}慢逾，常稱之為“擬合優(yōu)度”和“修正的擬合優(yōu)度”，指回歸方程對(duì)樣本的擬合程度幾何躏吊，這里可以看到氛改，修正的擬合優(yōu) 度=0.9104，表示擬合程度良好比伏，這個(gè)值當(dāng)然是越高越好胜卤。
5. F-statistic
   F統(tǒng)計(jì)量，也成為F檢驗(yàn)赁项，常常用于判斷方程整體的顯著性檢驗(yàn)葛躏，其值越大越顯著；其P值為p-value: < 2.2e-16悠菜，顯然是<0.05的舰攒，可以認(rèn)為方程在P=0.05的水平上還是通過(guò)顯著性檢驗(yàn)的。

注意到：Coefficients: (1 not defined because of singularities)
由于潛在的依賴關(guān)系導(dǎo)致有1個(gè)特征對(duì)輸出的作用無(wú)法和其他特征分清楚悔醋，這種現(xiàn)象：混疊（aliasing）摩窃。
使用alias函數(shù)顯示從模型中去除的特征。

alias(cars_model1) 
cars_model2 <- lm(Price~.-Saturn,data=cars_train) #去掉Saturn行
summary(cars_model2)

發(fā)現(xiàn)有問(wèn)題的是Saturn（代表車輛是否是土星）芬骄。

1）殘差分析

殘差：模型對(duì)特定觀測(cè)數(shù)據(jù)產(chǎn)生的誤差猾愿，輸出的實(shí)際值和預(yù)測(cè)值之間的差異???? = ???? ? ????^。
通常會(huì)把殘差從小到大排序账阻，較為理想的是0中位數(shù)和較小的四分位數(shù)值蒂秘。

可以看到二手車的均價(jià)在21k，但是50%的數(shù)在±1.6k左右淘太，這個(gè)結(jié)果較為合理姻僧。

殘差圖
R語(yǔ)言里的模型診斷圖（Residuals vs Fitted规丽，Normal QQ ， Scale-Location 撇贺，Residuals Leverage）
?分位圖（Quantile-Quantile plot, Q-Q plot）：
通過(guò)比較分位數(shù)（quantile）的值來(lái)比較兩種分布赌莺。

在理想線性模型中有五大假設(shè)。其中之一便是殘差應(yīng)該是一個(gè)正態(tài)分布显熏，與估計(jì)值無(wú)關(guān)雄嚣。如果殘差還和估計(jì)值有關(guān)，那就說(shuō)明模型仍然有值得去改進(jìn)的地方喘蟆，當(dāng)然缓升，這樣的模型準(zhǔn)確程度也就大打折扣。

QQ-plot用來(lái)檢測(cè)其殘差是否是正態(tài)分布的蕴轨。對(duì)應(yīng)于正態(tài)分布的QQ圖,就是由標(biāo)準(zhǔn)正態(tài)分布的分位數(shù)為橫坐標(biāo),樣本值為縱坐標(biāo)的散點(diǎn)圖港谊；橫坐標(biāo)也可以是1，2橙弱，3歧寺，表示幾個(gè)標(biāo)準(zhǔn)差。

# 殘差
machine_residuals <- machine_model1$residuals 
# lm線性方程擬合值
machine_fitted_values <- machine_model1$fitted.values 
# 179個(gè)訓(xùn)練集的序列號(hào)
machine_train_ids <- rownames(machine_train)
#離群殘差點(diǎn):大于150的,放入序列號(hào),否則為空;abs()取絕對(duì)值
machine_large_residuals <- ifelse(abs(machine_residuals) > 150,machine_train_ids,'')
 
p1 <- qplot(machine_fitted_values,machine_residuals) #橫坐標(biāo)是擬合值棘脐，縱坐標(biāo)是殘差
p1 <- p1 + ggtitle("Residual Plot for CPU Data Set") #標(biāo)題
p1 <- p1 + theme(plot.title = element_text(lineheight=.8, face="bold")) #字體
p1 <- p1 + xlab("Fitted Values")  #橫坐標(biāo)名稱
p1 <- p1 + ylab("Residuals") #縱坐標(biāo)名稱
p1 <- p1 + geom_text(size = 4, hjust=-0.15, vjust=0.1, aes(label=machine_large_residuals)) # 標(biāo)記離群殘差點(diǎn)
p1

cars_residuals <- cars_model1$residuals
cars_fitted_values <- cars_model1$fitted.values
cars_train_ids <- rownames(cars_train)
cars_large_residuals <- ifelse(abs(cars_residuals) > 9500,cars_train_ids,'')
 
p2 <- qplot(cars_fitted_values,cars_residuals) 
p2 <- p2 + ggtitle("Residual Plot for Cars Data Set") 
p2 <- p2 + theme(plot.title = element_text(lineheight=.8, face="bold")) 
p2 <- p2 + xlab("Fitted Values")  
p2 <- p2 + ylab("Residuals")
p2 <- p2 + geom_text(size = 4, hjust=-0.15, vjust=0.1, aes(label=cars_large_residuals))
p2

#par(mfrow=c(2,1)) #畫布布局斜筐，兩行一列
par(mar=c(1,1,1,1)) #設(shè)置圖形的邊界，下蛀缝，左顷链，上，右
# qqnorm生成y中值的正常QQ圖
qqnorm(machine_residuals, main = "Normal Q-Q Plot for CPU data set")
# qqline函數(shù)用于繪制QQ圖的近似擬合直線屈梁，其解析式a是正態(tài)分布的標(biāo)準(zhǔn)差嗤练，截距b為均值
qqline(machine_residuals)
qqnorm(cars_residuals, main = "Normal Q-Q Plot for Cars data set")
qqline(cars_residuals)

2）線性回歸的性能衡量指標(biāo)

?F統(tǒng)計(jì)量：來(lái)源于檢查兩個(gè)（正態(tài)）分布的方差之間是否存在統(tǒng)計(jì)顯著性的F檢驗(yàn)。
檢驗(yàn)只有截距的模型的殘差方差和現(xiàn)有訓(xùn)練模型的殘差方差的顯著性差異在讶。
anova()函數(shù)代表方差分析（analysis of variance）

# 概率分布的T分布煞抬，調(diào)用t分布P值函數(shù)pt()即可獲得該統(tǒng)計(jì)量的P值。
# 參數(shù)：q:t統(tǒng)計(jì)量的值,df:自由度,lower.tail:確定計(jì)算概率的方向构哺。如果lower.tail=T革答，計(jì)算Pr(X ≤x),反之,計(jì)算Pr(X >x)。
(q <- 5.210e-02 / 1.885e-02)
pt(q, df = 172, lower.tail = F) * 2
# 只有截距的模型
machine_model_null = lm(PRP~1,data=machine_train)
anova(machine_model_null, machine_model1)

n_machine <- nrow(machine_train)
k_machine <- length(machine_model1$coefficients) -1 #6列變量
# 殘差標(biāo)準(zhǔn)差RSE
sqrt(sum(machine_model1$residuals ^ 2) / (n_machine - k_machine - 1)) 
n_cars <- nrow(cars_train)
k_cars <- length(cars_model1$coefficients) -1
sqrt(sum(cars_model1$residuals ^ 2) / (n_cars - k_cars - 1))

mean(machine_train$PRP)
mean(cars_train$Price)

# 計(jì)算R^2
compute_rsquared <- function(x,y) {
  rss <- sum((x-y)^2)
  tss <- sum((y-mean(y))^2)
  return(1-(rss/tss))
}

compute_rsquared(machine_model1$fitted.values,machine_train$PRP)
compute_rsquared(cars_model2$fitted.values,cars_train$Price)

3）比較不同的回歸模型

?如何比較兩個(gè)輸入特征數(shù)量不相同的回歸模型曙强？
一個(gè)回歸模型是??1個(gè)特征數(shù)量残拐，另一個(gè)回歸模型是??2個(gè)特征數(shù)量
兩個(gè)模型該如何比較??^2？
調(diào)整后的adjusted ??^2

compute_adjusted_rsquared <- function(x,y,p) {
  n <- length(y)
  r2 <- compute_rsquared(x,y)
  return(1 - ((1 - r2) * (n-1)/(n-p-1)))
}

compute_adjusted_rsquared(machine_model1$fitted.values,machine_train$PRP,k_machine)
compute_adjusted_rsquared(cars_model2$fitted.values,cars_train$Price,k_cars)

4）在測(cè)試集上的性能

?在訓(xùn)練集和測(cè)試集都需要去觀察模型的性能
用訓(xùn)練集和測(cè)試集去計(jì)算性能（均方差）
函數(shù)predict

#創(chuàng)建一個(gè)函數(shù)旗扑，來(lái)計(jì)算通過(guò)模型得到的結(jié)果和實(shí)際結(jié)果之間的均方差
compute_mse <- function(predictions, actual) { 
  mean((predictions-actual)^2) }

machine_model1_predictions <- predict(machine_model1, machine_test)
compute_mse(machine_model1$fitted.values, machine_train$PRP)
compute_mse(machine_model1_predictions, machine_test$PRP)

cars_model2_predictions <- predict(cars_model2, cars_test)
compute_mse(cars_model2$fitted.values, cars_train$Price)
compute_mse(cars_model2_predictions, cars_test$Price)

5.線性回歸的問(wèn)題

1)多重共線性

是指線性回歸模型中解釋變量之間由于存在精確相關(guān)關(guān)系或高度相關(guān)關(guān)系而導(dǎo)致模型估計(jì)失真或難以估計(jì)蹦骑。
多重共線性如何觀測(cè)到慈省？
?可疑處一：
兩個(gè)高度共線性的特征具有較大的p值臀防，但是如果去掉其中一個(gè)眠菇，另一個(gè)p值變小
?可疑處二：
某個(gè)系數(shù)出現(xiàn)不正常的符號(hào)，如：預(yù)測(cè)收入的模型袱衷，教育背景的系數(shù)是負(fù)數(shù)
處理辦法
?把兩個(gè)特征合并
?直接去除其中一個(gè)

多重共線性可以對(duì)線性模型中的每個(gè)輸入特征計(jì)算其方差膨脹因子（variance inflation factor, VIF）捎废。
VIF的計(jì)算步驟為：
1.用某特征作為輸出，其他作為輸入計(jì)算??2致燥；

2.

在R中可直接使用VIF函數(shù)：

#在R中可直接使用VIF函數(shù)
library(carData)
library("car")
vif(cars_model2)

#計(jì)算sedan的VIF,用擬合值和真實(shí)值作為輸入變量
sedan_model <- lm(sedan ~.-Price-Saturn, data=cars_train)
sedan_r2 <- compute_rsquared(sedan_model$fitted.values,cars_train$sedan)
1 - (1-sedan_r2)

如果vif的值超過(guò)4或者更大的特征就是可疑的登疗；如果vif的值大于10就有多重共線性的極大可能性。

2）離群值

離群值可能由于測(cè)量誤差產(chǎn)生嫌蚤，沒(méi)有選對(duì)特征或創(chuàng)建了錯(cuò)誤的種類導(dǎo)致辐益。
?小心去除離群值
包含它們可能會(huì)產(chǎn)生顯著改變預(yù)測(cè)模型系數(shù)的效果

離群值一般通過(guò)殘差圖來(lái)看。殘差就是預(yù)測(cè)值和真實(shí)值之間的差異脱吱。殘差圖可以直接用plot畫出智政。

#離群值一般通過(guò)殘差圖來(lái)看。
plot(cars_model2)

#第200行的PRP的值很大箱蝠，去掉這個(gè)離群點(diǎn)
machine_model2 <- lm(PRP~.,
      data=machine_train[!(rownames(machine_train)) %in% c(200),])
summary(machine_model2)

machine_model2_predictions <- predict(machine_model2, machine_test)
compute_mse(machine_model2_predictions, machine_test$PRP)

說(shuō)明：去掉離群點(diǎn)续捂，model2計(jì)算出的MSE應(yīng)該比model1的小。

3)特征選擇

實(shí)際環(huán)境中特征數(shù)量太多宦搬，在模型中選擇一個(gè)特征子集以構(gòu)成一個(gè)帶有更少特征的新模型的過(guò)程牙瓢。
前向選擇（forward selection）
逐步回歸（step regression），從一個(gè)沒(méi)有特征的空模型開始间校，接著進(jìn)行??次（針對(duì)每個(gè)特征一次）簡(jiǎn)單線性回歸并從中選最好的矾克。
后向選擇（backward selection）
AIC指標(biāo)：越小越好

AIC準(zhǔn)則是由日本統(tǒng)計(jì)學(xué)家Akaike與1973年提出的，全稱是最小化信息量
準(zhǔn)則（Akaike Information Criterion）撇簿。它是擬合精度和參數(shù)個(gè)數(shù)的加權(quán)函數(shù)：
AIC=2（模型參數(shù)的個(gè)數(shù)）-2ln（模型的極大似然函數(shù)）

#建立一個(gè)只有截距項(xiàng)的模型
machine_model_null = lm(PRP~1,data=machine_train)
machine_model1 <- lm(PRP~.,data=machine_train)
#采用向前選擇方式來(lái)進(jìn)行特征選擇
machine_model3 <- step(machine_model_null, 
        scope = list(lower = machine_model_null, 
                     upper=machine_model1), 
        direction = "forward")

cars_model_null <- lm(Price~1,data=cars_train)
cars_model2 <- lm(Price~.-Saturn,data=cars_train)
#采用向后選擇方式來(lái)進(jìn)行特征選擇
cars_model3 <- step(cars_model2, 
                    scope=list(lower=cars_model_null, 
                               upper=cars_model2), 
                    direction="backward")

#預(yù)測(cè)
machine_model3_predictions <- predict(machine_model3, machine_test)
#評(píng)價(jià)模型
compute_mse(machine_model3_predictions, machine_test$PRP)

cars_model3_predictions <- predict(cars_model3, cars_test)
compute_mse(cars_model3_predictions, cars_test$Price)

6.正則化與過(guò)擬合

過(guò)擬合
模型在理想情況下聂渊，應(yīng)該是對(duì)訓(xùn)練集和測(cè)試集都能擬合得比較好。如果模型只對(duì)訓(xùn)練集擬合得很好四瘫，對(duì)測(cè)試集擬合得一般汉嗽，那就說(shuō)明存在過(guò)擬合得情況。模型的魯棒性較差找蜜，無(wú)法適應(yīng)普片情況饼暑。
對(duì)過(guò)擬合和欠擬合進(jìn)行直觀的理解：

嶺回歸（ridge regression）：通過(guò)其約束條件引入偏誤但能有效第減少模型的方差包雀。

所以，正則化的關(guān)鍵就是尋找合適的??房资。一般采用的是交叉驗(yàn)證法來(lái)確定??蜕劝。

在R中，可以使用glmnet包來(lái)實(shí)現(xiàn)嶺回歸和Lasso回歸轰异。
glmnet是由斯坦福大學(xué)的統(tǒng)計(jì)學(xué)家們開發(fā)的一款R包岖沛，用于在傳統(tǒng)的廣義線性回歸模型的基礎(chǔ)上添加正則項(xiàng)，以有效解決過(guò)擬合的問(wèn)題搭独，支持線性回歸婴削，邏輯回歸，泊松回歸牙肝，cox回歸等多種回歸模型唉俗。
glmet接受一個(gè)矩陣,每一行為一個(gè)觀測(cè)向量，每一列代表一個(gè)特征配椭。y是響應(yīng)變量虫溜。
alpha=0代表嶺回歸，alpha=1代表lasso回歸
alpah如果介于0和1之間股缸，則代表既有嶺回歸又與lasso回歸的混合模型——彈性網(wǎng)絡(luò)衡楞，此時(shí)aplha代表混合比。
對(duì)于每種模型Glmnet都提供了glmnet用于擬合模型, cv.glmnet使用k折交叉驗(yàn)證擬合模型, predict對(duì)數(shù)據(jù)進(jìn)行預(yù)測(cè)（分類/回歸）敦姻，coef用于提取指定lambda時(shí)特征的系數(shù)瘾境。

library(Matrix)
library(glmnet)
#使用model.matrix先來(lái)創(chuàng)建特征矩陣坎背，同時(shí)確保各列都是數(shù)值型（邏輯型、數(shù)值型寄雀、因子等）
cars_train_mat <- model.matrix(Price~.-Saturn, cars_train)[,-1]
#給定lamda范圍,從10^8 到 10^-4,平均生成250個(gè)lamda
lambdas <- 10 ^ seq(8,-4,length=250)

#嶺回歸
cars_models_ridge= glmnet(cars_train_mat,cars_train$Price,
                          alpha=0,lambda=lambdas)
#lasso回歸
cars_models_lasso= glmnet(cars_train_mat,cars_train$Price,
                          alpha=1,lambda=lambdas)

cars_models_ridge
#選出第70個(gè)lambda
cars_models_ridge$lambda[70]
#選出第70個(gè)模型的系數(shù)
coef(cars_models_ridge)[,70]

#畫出250個(gè)模型的系數(shù)隨著lamda的變化而變化。橫坐標(biāo)是lamda陨献，縱坐標(biāo)是各特征參數(shù)的取值
layout(matrix(c(1,2), 1, 2))
plot(cars_models_ridge, xvar = "lambda", main = "Coefficient Values vs. Log Lambda for Ridge Regression")
plot(cars_models_lasso, xvar = "lambda", main = "Coefficient Values vs. Log Lambda for Lasso")

可看到盒犹，隨著??的增大，各特征的系數(shù)被壓縮到了0眨业。此時(shí)會(huì)導(dǎo)致欠擬合急膀。那么??選多大合適呢？可以使用交叉檢驗(yàn)法來(lái)選擇合適的??龄捡。
R中卓嫂，可以直接使用cv.glmnet 來(lái)幫忙選擇最優(yōu)的??。

ridge.cv <- cv.glmnet(cars_train_mat,cars_train$Price,
                      alpha=0,lambda=lambdas)
lambda_ridge <- ridge.cv$lambda.min
lambda_ridge

lasso.cv <- cv.glmnet(cars_train_mat,cars_train$Price,
                      alpha=1,lambda=lambdas)
lambda_lasso <- lasso.cv$lambda.min
lambda_lasso

#x軸代表經(jīng)過(guò)log以后的lambda值聘殖，y軸代表模型的誤差晨雳，cv.glmnet會(huì)自動(dòng)選擇使誤差最小的lambda（左側(cè)的虛線）
layout(matrix(c(1,2), 1, 2))
plot(ridge.cv, col = gray.colors(1))
title("Ridge Regression", line = +2)
plot(lasso.cv, col = gray.colors(1))
title("Lasso", line = +2)

同時(shí)也可以使用coef提取每一個(gè)特征在指定lambda下的系數(shù)：

#提取lambda = 9.201432時(shí)的特征系數(shù)，·代表經(jīng)過(guò)L1正則化后這些特征已經(jīng)被消掉了奸腺。
coef.apprx = coef(ridge.cv, s = 9.201432)
coef.apprx

有了??后餐禁，就可以來(lái)進(jìn)行預(yù)測(cè)。

#輸出新數(shù)據(jù)的預(yù)測(cè)值突照，type參數(shù)允許選擇預(yù)測(cè)的類型并提供預(yù)測(cè)值帮非，newx代表要預(yù)測(cè)的數(shù)據(jù)。
predict(cars_models_lasso, type="coefficients", s = lambda_lasso)

cars_test_mat <- model.matrix(Price~.-Saturn, cars_test)[,-1]
cars_ridge_predictions <- predict(cars_models_ridge, s = lambda_ridge, newx = cars_test_mat)
compute_mse(cars_ridge_predictions, cars_test$Price)
cars_lasso_predictions <- predict(cars_models_lasso, s = lambda_lasso, newx = cars_test_mat)
compute_mse(cars_lasso_predictions, cars_test$Price)

【完整步驟】從機(jī)器學(xué)習(xí)的角度來(lái)做線性回歸：

第一步：數(shù)據(jù)預(yù)處理：

#二手車交易數(shù)據(jù)
data(cars)
#進(jìn)行數(shù)據(jù)預(yù)處理讹蘑，查看一下各特征之間的相關(guān)系數(shù)矩陣末盔，判斷相關(guān)關(guān)系
cars_cor <- cor(cars[-1]) #去掉第一列后，進(jìn)行相關(guān)系數(shù)統(tǒng)計(jì)
findCorrelation(cars_cor)
findCorrelation(cars_cor, cutoff = 0.75)
#通過(guò)查看相關(guān)矩陣座慰，發(fā)現(xiàn)在Doors和coupe之間存在較高的相關(guān)性：如果是coupe很有可能是2雙門陨舱，否則是4門
cor(cars$Doors,cars$coupe)
table(cars$coupe,cars$Doors)  #用交叉列聯(lián)表來(lái)查看相關(guān)性
#查找完全線性組合，發(fā)現(xiàn)15列和18列存在完全線性組合
findLinearCombos(cars)  
#根據(jù)建議版仔，去掉具有完全線性組合的特征
cars <- cars[,c(-15,-18)]

第二步：劃分訓(xùn)練集和測(cè)試集

#劃分訓(xùn)練集和測(cè)試集隅忿，同時(shí)標(biāo)注特征數(shù)據(jù)和標(biāo)簽數(shù)據(jù)
cars_sampling_vector <- createDataPartition(cars$Price, p=0.85, list = FALSE)

cars_train <- cars[cars_sampling_vector,]
cars_train_features <- cars[,-1]
cars_train_labels <- cars$Price[cars_sampling_vector] 

cars_test <- cars[-cars_sampling_vector,]
cars_test_labels <- cars$Price[-cars_sampling_vector] 

第三步：訓(xùn)練模型：

#根據(jù)訓(xùn)練集建立模型
cars_model1 <- lm(Price~.,data=cars_train)

第四步：模型優(yōu)化
模型評(píng)估（包括：殘差分析、顯著性檢驗(yàn)等邦尊。其中顯著性檢驗(yàn)又包括線性關(guān)系檢驗(yàn)和回歸系數(shù)檢驗(yàn)）

summary(cars_model1)
#混疊(aliasing),Coefficients(1 not defined because of singularities),Saturn行數(shù)據(jù)都為NA
#移除特征再進(jìn)行回歸
cars_model2 <- lm(Price~.-Saturn,data=cars_train)
summary(cars_model2)
#殘差分析
cars_residuals <- cars_model2$residuals
qqnorm(cars_residuals, main = "Normal Q-Q Plot for Cars data set")
qqline(cars_residuals)
#顯著性檢驗(yàn)
#方差分析
cars_model_null = lm(Price~1,data=cars_train)
anova(cars_model_null, cars_model2)
#計(jì)算R^2
compute_rsquared <- function(x,y) {
  rss <- sum((x-y)^2)
  tss <- sum((y-mean(y))^2)
  return(1-(rss/tss))
}
compute_rsquared(cars_model2$fitted.values,cars_train$Price)
#調(diào)整后的adjusted ??^2
compute_adjusted_rsquared <- function(x,y,p) {
  n <- length(y)
  r2 <- compute_rsquared(x,y)
  return(1 - ((1 - r2) * (n-1)/(n-p-1)))
}
k_cars <- length(cars_model2$coefficients) -1
compute_adjusted_rsquared(cars_model2$fitted.values,cars_train$Price,k_cars)

第五步：預(yù)測(cè)

cars_model2_predictions <- predict(cars_model2, cars_test)

第六步：模型評(píng)價(jià)
通過(guò)分別計(jì)算訓(xùn)練集上的模型和測(cè)試集上的模型得到結(jié)果與實(shí)際結(jié)果的差值比較來(lái)判斷模型的好壞（MSE）背桐。
越小的MSE或者RSS，代表擬合程度越好蝉揍。

#創(chuàng)建一個(gè)函數(shù)链峭，來(lái)計(jì)算通過(guò)模型得到的結(jié)果和實(shí)際結(jié)果直接的均方差
compute_mse <- function(predictions, actual) { 
  mean((predictions-actual)^2) 
}
compute_mse(cars_model2$fitted.values, cars_train$Price)
compute_mse(cars_model2_predictions, cars_test$Price)