1.线性回归入门

线性回归的假设

假设输出变量是一组特征变量的加权线性函数

◾假设1：误差项具有同方差性

同方差（homoscedasticity）：误差部分的方差不会随着输入特征的值或水平而变化
异方差（heteroskedastic）：误差部分的方差随着输入特征的值的变化而变化

order()函数用于返回向量大小顺序的秩，从小到大的数字的序号值。
x<-c(10,6,4,7,8)
order(x)
[1] 3 2 4 5 1
sapply(x,FUN)，第一个参数x是需要处理的数据，FUN是处理数据x的Funtion。

library(ggplot2)
set.seed(15)
# 随机产生200个数，均值100，标准差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)})
#标准差随x变大而变大，y的值也随之变大，误差项的方差与x有相关性。其实是构造了一个随x变化的异方差的y的数据
df <- data.frame(x1 = x1, y = y)
fit<-coef(lm(y ~ x1, data = df))
# coef()函数提取线性回归lm模型的系数
# qplot画点图，geom_abline画线图
# size：坐标轴标签字体大小，lineheight: 标签行间距的倍数，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

可以看到，系数是17.0848，1.81745，非常接近我们构造的函数17+1.8*x + rnorm(1,mean=0, sd = x/4)，因为有噪声，所以截距和x的系数有细微变化。

同方差性与异方差性的区别：

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')

◾假设2：误差项具有独立性

人工扩大数据集的方式会破坏独立性的假设。
例如：取出一部分数据，把特征值和输出值都乘以2，看起来数据集合变大了，但是数据之间的误差有了相关性。
不同工厂的测量数据，工厂应该作为一种特征。

◾假设3：特征在统计学上是相互独立的

2.简单线性回归

从均匀分布中随机抽取点，让数据理想地分开
𝑦 = 1.67𝑥1 − 2.93 + 𝑁(0,2^2)

估算回归系数

注意：
如果 𝛽1^ = 0表示输出变量与输入特征相互独立
即使两个变量统计学上相互独立，通常还是会有少量的协方差
协方差是0也不代表是相互独立
如果训练一个线性回归模型，一般𝛽1^ ≠ 0

set.seed(5427395)
nObs <- 100
x1minrange <- 5 
x1maxrange <- 25
# 随机产生100个最小值为5，最大值为25的均匀分布随机数
x1 <- runif(nObs,x1minrange,x1maxrange)
# 随机产生100个均值为0，标准差为2的正态分布随机数
e <- rnorm(nObs,mean = 0, sd = 2.0)
y <- 1.67*x1 - 2.93 + e
# 构造一个数据框格式
df <- data.frame(y,x1)
# 线性拟合
myfit <- lm(y~x1,df)
summary(myfit)

因为e的噪声干扰，发现截距在-2.93附近，x1自变量系数在1.67附近。e是同方差干扰。

# 画图
p <- qplot(x1,y,data=df) 
p <- p + ggtitle("Simple Linear Regression")  #画标题
# 画lm函数拟合直线
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固定系数下原始直线
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

以上是回归线性系数估计拟合线与固定系数和截距划线的比较，发现两条线基本重合，拟合度较好。

3.多元线性回归

具有k个输入特征的多元线性回归模型𝑦 = 𝛽𝑘𝑥𝑘 + 𝛽𝑘−1𝑥𝑘−1 + ⋯ + 𝛽1𝑥1 + 𝛽0 + 𝜀
通常假设𝜀与(𝑥1, 𝑥2, … , 𝑥𝑘)相互独立。

1）预测CPU性能

把数据集做成训练集和测试集，比例85%:15%
1.该数据没有缺失值，而且样本较少，仅有200行

library(caret)
set.seed(4352345)
machine$PRP
# 按照输出变量PRP进行索引，提取85%作为训练集
machine_sampling_vector <- createDataPartition(machine$PRP, p = 0.85, list = FALSE)
# 获取训练集数据
machine_train <- machine[machine_sampling_vector,]
machine_train_labels <- machine$PRP[machine_sampling_vector]

# 获取测试集数据
machine_test <- machine[-machine_sampling_vector,]
machine_test_labels <- machine$PRP[-machine_sampling_vector]

2.数据设置好后，分析和检查针对线性回归的假设

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

# 筛选与其他相关性系数强的变量并且删除
findCorrelation(machine_correlations) 
findCorrelation(machine_correlations, cutoff = 0.75)
cor(machine_train$MMIN,machine_train$MMAX)

高度相关性的阈值默认的0.9没有发现哪个特征该去掉。
但是，当阈值改成0.75的时候，caret包推荐去掉第三个特征MMAX。
检查发现MMIN和MMAX相关性很大。最小主内存的型号也有较大的主内存，这里我们不修改。

2）预测二手汽车的价格

数据集caret，获取二手车可靠价格的数据，包含804种通用汽车（GM）品牌的汽车。

第一步：数据预处理

library(caret)
data(cars)
# 查看一下各特征之间的相关系数矩阵，判断相关关系
cars_cor <- cor(cars[,-1]) #去掉第一列price
findCorrelation(cars_cor) #默认cutoff=0.9
findCorrelation(cars_cor, cutoff=0.75)
cor(cars$Doors,cars$coupe)
# 用交叉列联表来查看相关性
table(cars$coupe,cars$Doors)
# 查找完全线性组合的特征
findLinearCombos(cars)
# 根据建议，去掉具有完全线性组合的特征
cars <- cars[,c(-15,-18)] 

第三列doors有超过75%相关性。查看相关矩阵，发现在Doors和coupe之间存在较高的相关性。
跑车大概率是2门车，跑车和车门相关性为-0.8254435。门越少，是跑车可能性越高。2门车中140个是跑车，50个不是跑车，4门车中0个是跑车。

函数结果表示：要把15(coupe)和18(wagon)remove掉，因为15和18是其他特征的组合。

第二步：把数据分为训练集和测试集

set.seed(232455)
# 划分训练集和测试集，同时标注特征数据和标签数据
cars_sampling_vector <- createDataPartition(cars$Price, p = 0.85, list = FALSE)
# 获取训练集数据
cars_train <- cars[cars_sampling_vector,]
cars_train_features <- cars[,-1] #去掉第一列price
cars_train_labels <- cars$Price[cars_sampling_vector]
# 获取测试集数据
cars_test <- cars[-cars_sampling_vector,]
cars_test_labels <- cars$Price[-cars_sampling_vector] 

4.评估线性回归模型

利用lm函数，用线性回归模型来拟合数据。

# 根据训练集建立模型
machine_model1 <- lm(PRP~.,data=machine_train)
cars_model1 <- lm(Price~.,data=cars_train)
# 模型评估（包括：残差分析、显著性检验等。其中显著性检验又包括线性关系检验和回归系数检验）
summary(cars_model1)
summary(machine_model1)

summary()结果解读

1. 调用：Call，表明lm是如何被调用的
2. 残差统计量：Residuals，残差第一四分位数（1Q）和第三分位数（Q3）有大约相同的幅度，意味着有较对称的钟形分布
3. 系数：Coefficients
   分别表示： 估值标准误差 T值 P值
   Intercept：表示截距
   Mileage：影响因子/特征
   Estimate：包含由普通最小二乘法计算出来的估计回归系数。
   Std. Error：估计的回归系数的标准误差。
   t值：系数估计值与标准误差的比值。
   P值：估计系数不显著的可能性，有较大P值的变量是可以从模型中移除的候选变量。可以直接通过P值与预设的0.05进行比较，来判定对应的解释变量的显著性，检验的原假设是：该系数显著为0；若P<0.05，则拒绝原假设，即对应的变量显著不为0。（p值很高并不一定代表特征值和输出没有线性关系；它只表明有其他特征存在的时候，这个特征对输出变量不提供新的信息）
   可以看到Mileage、Cylinder、Doors都可以认为是在P为0.05的水平下显著不为0，通过显著性检验；Cruise的P值为0.20025，不显著。
4. Multiple R-squared和Adjusted R-squared
   这两个值，即R^{2}，常称之为“拟合优度”和“修正的拟合优度”，指回归方程对样本的拟合程度几何，这里可以看到，修正的拟合优 度=0.9104，表示拟合程度良好，这个值当然是越高越好。
5. F-statistic
   F统计量，也成为F检验，常常用于判断方程整体的显著性检验，其值越大越显著；其P值为p-value: < 2.2e-16，显然是<0.05的，可以认为方程在P=0.05的水平上还是通过显著性检验的。

注意到：Coefficients: (1 not defined because of singularities)
由于潜在的依赖关系导致有1个特征对输出的作用无法和其他特征分清楚，这种现象：混叠（aliasing）。
使用alias函数显示从模型中去除的特征。

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

发现有问题的是Saturn（代表车辆是否是土星）。

1）残差分析

残差：模型对特定观测数据产生的误差，输出的实际值和预测值之间的差异𝑒𝑖 = 𝑦𝑖 − 𝑦𝑖^。
通常会把残差从小到大排序，较为理想的是0中位数和较小的四分位数值。

可以看到二手车的均价在21k，但是50%的数在±1.6k左右，这个结果较为合理。

残差图
R语言里的模型诊断图（Residuals vs Fitted，Normal QQ ， Scale-Location ，Residuals Leverage）
◾分位图（Quantile-Quantile plot, Q-Q plot）：
通过比较分位数（quantile）的值来比较两种分布。

在理想线性模型中有五大假设。其中之一便是残差应该是一个正态分布，与估计值无关。如果残差还和估计值有关，那就说明模型仍然有值得去改进的地方，当然，这样的模型准确程度也就大打折扣。

QQ-plot用来检测其残差是否是正态分布的。对应于正态分布的QQ图,就是由标准正态分布的分位数为横坐标,样本值为纵坐标的散点图；横坐标也可以是1，2，3，表示几个标准差。

# 残差
machine_residuals <- machine_model1$residuals 
# lm线性方程拟合值
machine_fitted_values <- machine_model1$fitted.values 
# 179个训练集的序列号
machine_train_ids <- rownames(machine_train)
#离群残差点:大于150的,放入序列号,否则为空;abs()取绝对值
machine_large_residuals <- ifelse(abs(machine_residuals) > 150,machine_train_ids,'')
 
p1 <- qplot(machine_fitted_values,machine_residuals) #横坐标是拟合值，纵坐标是残差
p1 <- p1 + ggtitle("Residual Plot for CPU Data Set") #标题
p1 <- p1 + theme(plot.title = element_text(lineheight=.8, face="bold")) #字体
p1 <- p1 + xlab("Fitted Values")  #横坐标名称
p1 <- p1 + ylab("Residuals") #纵坐标名称
p1 <- p1 + geom_text(size = 4, hjust=-0.15, vjust=0.1, aes(label=machine_large_residuals)) # 标记离群残差点
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)) #设置图形的边界，下，左，上，右
# qqnorm生成y中值的正常QQ图
qqnorm(machine_residuals, main = "Normal Q-Q Plot for CPU data set")
# qqline函数用于绘制QQ图的近似拟合直线，其解析式a是正态分布的标准差，截距b为均值
qqline(machine_residuals)
qqnorm(cars_residuals, main = "Normal Q-Q Plot for Cars data set")
qqline(cars_residuals)

2）线性回归的性能衡量指标

◾F统计量：来源于检查两个（正态）分布的方差之间是否存在统计显著性的F检验。
检验只有截距的模型的残差方差和现有训练模型的残差方差的显著性差异。
anova()函数代表方差分析（analysis of variance）

# 概率分布的T分布，调用t分布P值函数pt()即可获得该统计量的P值。
# 参数：q:t统计量的值,df:自由度,lower.tail:确定计算概率的方向。如果lower.tail=T，计算Pr(X ≤x),反之,计算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列变量
# 残差标准差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)

# 计算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）比较不同的回归模型

◾如何比较两个输入特征数量不相同的回归模型？
一个回归模型是𝑘1个特征数量，另一个回归模型是𝑘2个特征数量
两个模型该如何比较𝑅^2？
调整后的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）在测试集上的性能

◾在训练集和测试集都需要去观察模型的性能
用训练集和测试集去计算性能（均方差）
函数predict

#创建一个函数，来计算通过模型得到的结果和实际结果之间的均方差
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.线性回归的问题

1)多重共线性

是指线性回归模型中解释变量之间由于存在精确相关关系或高度相关关系而导致模型估计失真或难以估计。
多重共线性如何观测到？
◾可疑处一：
两个高度共线性的特征具有较大的p值，但是如果去掉其中一个，另一个p值变小
◾可疑处二：
某个系数出现不正常的符号，如：预测收入的模型，教育背景的系数是负数
处理办法
◾把两个特征合并
◾直接去除其中一个

多重共线性可以对线性模型中的每个输入特征计算其方差膨胀因子（variance inflation factor, VIF）。
VIF的计算步骤为：
1.用某特征作为输出，其他作为输入计算𝑅2；

2.

在R中可直接使用VIF函数：

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

#计算sedan的VIF,用拟合值和真实值作为输入变量
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的值超过4或者更大的特征就是可疑的；如果vif的值大于10就有多重共线性的极大可能性。

2）离群值

离群值可能由于测量误差产生，没有选对特征或创建了错误的种类导致。
◾小心去除离群值
包含它们可能会产生显著改变预测模型系数的效果

离群值一般通过残差图来看。残差就是预测值和真实值之间的差异。残差图可以直接用plot画出。

#离群值一般通过残差图来看。
plot(cars_model2)

#第200行的PRP的值很大，去掉这个离群点
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)

说明：去掉离群点，model2计算出的MSE应该比model1的小。

3)特征选择

实际环境中特征数量太多，在模型中选择一个特征子集以构成一个带有更少特征的新模型的过程。
前向选择（forward selection）
逐步回归（step regression），从一个没有特征的空模型开始，接着进行𝑘次（针对每个特征一次）简单线性回归并从中选最好的。
后向选择（backward selection）
AIC指标：越小越好

AIC准则是由日本统计学家Akaike与1973年提出的，全称是最小化信息量
准则（Akaike Information Criterion）。它是拟合精度和参数个数的加权函数：
AIC=2（模型参数的个数）-2ln（模型的极大似然函数）

#建立一个只有截距项的模型
machine_model_null = lm(PRP~1,data=machine_train)
machine_model1 <- lm(PRP~.,data=machine_train)
#采用向前选择方式来进行特征选择
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)
#采用向后选择方式来进行特征选择
cars_model3 <- step(cars_model2, 
                    scope=list(lower=cars_model_null, 
                               upper=cars_model2), 
                    direction="backward")

#预测
machine_model3_predictions <- predict(machine_model3, machine_test)
#评价模型
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.正则化与过拟合

过拟合
模型在理想情况下，应该是对训练集和测试集都能拟合得比较好。如果模型只对训练集拟合得很好，对测试集拟合得一般，那就说明存在过拟合得情况。模型的鲁棒性较差，无法适应普片情况。
对过拟合和欠拟合进行直观的理解：

岭回归（ridge regression）：通过其约束条件引入偏误但能有效第减少模型的方差。

所以，正则化的关键就是寻找合适的𝜆。一般采用的是交叉验证法来确定𝜆。

在R中，可以使用glmnet包来实现岭回归和Lasso回归。
glmnet是由斯坦福大学的统计学家们开发的一款R包，用于在传统的广义线性回归模型的基础上添加正则项，以有效解决过拟合的问题，支持线性回归，逻辑回归，泊松回归，cox回归等多种回归模型。
glmet接受一个矩阵,每一行为一个观测向量，每一列代表一个特征。y是响应变量。
alpha=0代表岭回归，alpha=1代表lasso回归
alpah如果介于0和1之间，则代表既有岭回归又与lasso回归的混合模型——弹性网络，此时aplha代表混合比。
对于每种模型Glmnet都提供了glmnet用于拟合模型, cv.glmnet使用k折交叉验证拟合模型, predict对数据进行预测（分类/回归），coef用于提取指定lambda时特征的系数。

library(Matrix)
library(glmnet)
#使用model.matrix先来创建特征矩阵，同时确保各列都是数值型（逻辑型、数值型、因子等）
cars_train_mat <- model.matrix(Price~.-Saturn, cars_train)[,-1]
#给定lamda范围,从10^8 到 10^-4,平均生成250个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个lambda
cars_models_ridge$lambda[70]
#选出第70个模型的系数
coef(cars_models_ridge)[,70]

#画出250个模型的系数随着lamda的变化而变化。横坐标是lamda，纵坐标是各特征参数的取值
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")

可看到，随着𝜆的增大，各特征的系数被压缩到了0。此时会导致欠拟合。那么𝜆选多大合适呢？可以使用交叉检验法来选择合适的𝜆。
R中，可以直接使用cv.glmnet 来帮忙选择最优的𝜆。

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轴代表经过log以后的lambda值，y轴代表模型的误差，cv.glmnet会自动选择使误差最小的lambda（左侧的虚线）
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)

同时也可以使用coef提取每一个特征在指定lambda下的系数：

#提取lambda = 9.201432时的特征系数，·代表经过L1正则化后这些特征已经被消掉了。
coef.apprx = coef(ridge.cv, s = 9.201432)
coef.apprx

有了𝜆后，就可以来进行预测。

#输出新数据的预测值，type参数允许选择预测的类型并提供预测值，newx代表要预测的数据。
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)

【完整步骤】从机器学习的角度来做线性回归：

第一步：数据预处理：

#二手车交易数据
data(cars)
#进行数据预处理，查看一下各特征之间的相关系数矩阵，判断相关关系
cars_cor <- cor(cars[-1]) #去掉第一列后，进行相关系数统计
findCorrelation(cars_cor)
findCorrelation(cars_cor, cutoff = 0.75)
#通过查看相关矩阵，发现在Doors和coupe之间存在较高的相关性：如果是coupe很有可能是2双门，否则是4门
cor(cars$Doors,cars$coupe)
table(cars$coupe,cars$Doors)  #用交叉列联表来查看相关性
#查找完全线性组合，发现15列和18列存在完全线性组合
findLinearCombos(cars)  
#根据建议，去掉具有完全线性组合的特征
cars <- cars[,c(-15,-18)]

第二步：划分训练集和测试集

#划分训练集和测试集，同时标注特征数据和标签数据
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] 

第三步：训练模型：

#根据训练集建立模型
cars_model1 <- lm(Price~.,data=cars_train)

第四步：模型优化
模型评估（包括：残差分析、显著性检验等。其中显著性检验又包括线性关系检验和回归系数检验）

summary(cars_model1)
#混叠(aliasing),Coefficients(1 not defined because of singularities),Saturn行数据都为NA
#移除特征再进行回归
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)
#显著性检验
#方差分析
cars_model_null = lm(Price~1,data=cars_train)
anova(cars_model_null, cars_model2)
#计算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)
#调整后的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)

第五步：预测

cars_model2_predictions <- predict(cars_model2, cars_test)

第六步：模型评价
通过分别计算训练集上的模型和测试集上的模型得到结果与实际结果的差值比较来判断模型的好坏（MSE）。
越小的MSE或者RSS，代表拟合程度越好。

#创建一个函数，来计算通过模型得到的结果和实际结果直接的均方差
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)