Regression

回归分析

概念：

• 回归分析属于统计学的基本模型，回归分析(Regression Analysis)是用来确定2个或2个以上变量间关系的一种统计分析方法。
• 在回归分析中，变量有2类：因变量\(Y\)和自变量\(X\)。因变量通常是指实际问题中所关心的指标，用\(Y\)表示。而自变量是影响因变量取值的一个变量，用\(X\)表示，如果有多个自变量则表示为\(X_1, X_2, \ldots, X_p\)。

回归分析研究的主要步骤：

• 确定因变量\(Y\)与自变量\(X_1, X_2, \ldots, X_n\)之间的定量关系表达式，即回归方程;
• 对回归方程的置信度进行检验;
• 判断自变量\(X_n(n=1,2,\ldots,m)\)对因变量\(Y\)的影响;
• 利用回归方程进行预测.

线性回归模型常用的R函数：

函数	用途
lm()	拟合回归模型
summary()	展示拟合模型的详细结果
coefficients()	列出拟合模型的模型参数（截距项和斜率）
confint()	提供模型参数的置信区间（默认95%）
fitted()	列出拟合模型的预测值
residuals()	列出拟合模型的残差值
anova()	生成一个拟合模型的方差分析表，或者比较两个或更多拟合模型的方差分析表
vcov()	列出模型参数的协方差矩阵
AIC()	输出赤池信息统计量
plot()	生成评价拟合模型拟合情况的四幅诊断图（用于回归诊断）
predict()	用拟合模型对新的数据集预测相应变量

回归诊断R函数

基础安装

函数	用途
plot()	生成评价拟合模型拟合情况的四幅诊断图（用于回归诊断）

car package

函数	用途
qqplot()	分位数比较图（正态性检验）
durbinWatsonTest()	对误差自相关性做Durbin-Watson检验（独立性检验）
crPlots()	成分与残差图（线性检验）
ncvTest()	对非恒定的误差方差做得分检验（异方差性检验）
spreadLevelPlot()	分散水平检验（异方差性检验）
outlierTest()	Bonferroni离群点检验（离群点检验）
avPlots()	添加的变量图形
inluencePlot()	回归影响图
scatterplot()	增强的散点图
scatterplotMatrix()	增强的散点图矩阵
vif()	方差膨胀因子（多重共线性检验，如果sqrt(VIF) > 2,则说明存在多重共线性）

gvlma package

函数	用途
gvlma()	对线性模型提供了一个单独的综合检验（通过/不通过）

一个生成学生化残差柱状图，并添加正态曲线、核密度曲线和轴须图的函数

residplot = function(fit, nbreaks = 10) {
    
    z = rstudent(fit)
    hist(z, breaks = nbreaks, freq = FALSE,
         xlab = "Studentized Residual", main = "Distribution of Errors")
    rug(jitter(z), col = "brown")
    curve(dnorm(x, mean = mean(z), sd = sd(z)), add = TRUE, col = "blue", lwd = 2)
    lines(density(z)$x, density(z)$y, col = "red", lwd = 2, lty = 2)
    legend("topright", legend = c("Normal Curve", "Kernel Density Curve"),
           lty = 1:2, col = c("blue", "red"), cex = 0.7)
    
}

多项式回归

函数	用途
lm()	拟合回归模型
I()	从算术的角度来解释括号中的元素，可用于多项式回归中
poly()	多项式回归
scatterplot()(car pachage)	功能强劲的的散点图函数：散点图、线性拟合曲线、平滑(loess)曲线、箱线图等

多元线性回归

函数	用途
lm()	拟合回归模型
I()	从算术的角度来解释括号中的元素，可用于多项式回归中
cor()	检查变量之间的相关性
scatterplotMatrix()	功能强劲的的散点图矩阵函数：散点图、线性拟合曲线、平滑(loess)曲线、密度图和轴须图等(car package)
effect()	交互项展示(effect package)

模型选择

模型比较

考虑模型预测精度（模型尽可能地拟合数据）和模型简洁度（一个简单且能复制的模型）的调和

函数	用途
anova()	比较两个嵌套模型的拟合优度
AIC()	考虑了统计拟合度以及拟合的参数数目

变量选择

• 逐步回归(stepwise method)
• 向前逐步回归(forward stepwise)
• 向后逐步回归(backward stepwise)
• 向前向后逐步回归(stepwise stepwise)
• 全子集回归(all-subsets regression)
• 岭回归(ridge regression)
• Lasso和最小角回归()
• 弹性网()
• 主成分回归()
• 偏最小二乘回归()

函数	用途
step()	实现逐步回归模型的方法，依据的是精确的AIC准则
stepAIC()	实现逐步回归模型的方法，依据的是精确的AIC准则(MASS package)
regsubsets()	通过R方、调整R方,调整的R方或Mallows Cp统计量等准则运用全子集回归来选择“最佳”模型(leaps package)
plot()	绘制全子集回归的结果(leaps package)
subsets()	绘制全子集回归的结果(car package)
lm.ridg()	岭回归(MASS package)
lars()	lasso和最小角回归(lars package)
glmnet()	lasso(glmnet package)
enet()	弹性网络(elasticnet package)
pcr()	主成分回归(pls package)
plsr()	偏最小二乘回归(pls package)

评价模型泛化能力

交叉验证

函数	用途
crossval()	实现K折交叉验证(bootstrap package)
shrinkage()	R方的交叉验证(自变函数)

shrinkage function

shrinkage = function(fit, k = 10) {
    
    require(bootstrap)
    theta.fit = function(x, y) {lsfit(x, y)} # Find the Least Squares Fit
    theta.predict = function(fit, x) {cbind(1, x) %*% fit$coef}
    x = fit$model[, 2:ncol(fit$model)]
    y = fit$model[, 1]
    results = crossval(x, y, theta.fit, theta.predict, ngroup = k)
    r2 = cor(y, fit$fitted.values) ^ 2
    r2cv = cor(y, results$cv.fit) ^ 2
    
    cat("Original R-square =", r2, "\n")
    cat(k, "Fold Cross-Validatef R-square =", r2cv, "\n")
    cat("Change =", r2 - r2cv, "\n")
    
}

变量相对重要性

变量重要性排序

函数	用途
cor()	计算预测变量与响应变量的相关系数（预测变量间不相关）
as.data.frame(scale(data))	比较标准化的回归系数
relaimpo package	dsadadada
relweights()	相对权重：对所有可能子模型添加一个预测变量引起的R平方增加量的一个近似值（自变函数）

relweights function

relweights = function(fit,...) {
    R = cor(fit$model); nvar = ncol(R)
    rxx = R[2:nvar, 2:nvar]; rxy = R[2:nvar, 1]
    svd = eigen(rxx); evec = svd$vectors
    ev = svd$values; delta = diag(sqrt(ev))
    lambda = evec %*% delta %*% t(evec); lambdasq = lambda ^ 2
    beta = solve(lambda) %*% rxy; rsquare = colSums(beta ^ 2)
    rawwgt = lambdasq %*% beta ^ 2; import = (rawwgt / rsquare) * 100
    lbls = names(fit$model[2:nvar])
    rownames(import) = lbls; colnames(import) = "Weights"
    barplot(t(import), names.arg = lbls,
            ylab = "% of R-Square",
            xlab = "Predictor Variables", 
            main = "Relative Importance of Predictor Variables", 
            sub = paste("R-Square=", round(rsquare, digits = 3)), ...)
    return(import)
}

1. 简单线性回归(一元线性回归)

一元线性回归分析是处理两个变量之间关系的最简单模型，是两个变量之间的线性相关关系。

1.1 目录:

• 一元线性回归介绍
• 数据集介绍
• 建立回归模型
• 回归参数估计
• 回归方程的显著性检验
• 残差分析和异常点检测
• 模型预测

1.1 一元线性回归介绍

概念：

• 如果回归分析中，只包括一个自变量\(X\)和一个因变量\(Y\)时，且它们的关系是线性的，那么这种回归分析称为一元线性回归分析。

1.2 数据集介绍

数据集为2016年3月1日,白天开盘的交易数据,为锌的2个期货合约的分钟线的价格数据。数据集包括有3列，索引列为时间，\(zn1.Close\)为ZN1604合约的1分钟线的报价数据，\(zn2.Close\)为ZN1605合约的1分钟线的报价数据。

library(readr)
sim_linear <- read_csv("E:/R/Data Mining/Regression/data_image/simple_linear_regression.txt",
                         col_names = c("date", "x", "y"))

head(sim_linear);str(sim_linear)

##                  date     x     y
## 1 2016-03-01 09:01:00 14075 14145
## 2 2016-03-01 09:02:00 14095 14160
## 3 2016-03-01 09:03:00 14095 14160
## 4 2016-03-01 09:04:00 14095 14165
## 5 2016-03-01 09:05:00 14120 14190
## 6 2016-03-01 09:06:00 14115 14180

## Classes 'tbl_df', 'tbl' and 'data.frame':    27 obs. of  3 variables:
##  $ date: POSIXct, format: "2016-03-01 09:01:00" "2016-03-01 09:02:00" ...
##  $ x   : int  14075 14095 14095 14095 14120 14115 14110 14110 14105 14105 ...
##  $ y   : int  14145 14160 14160 14165 14190 14180 14170 14175 14170 14170 ...

利用散点图探索变量间的关系：

library(ggplot2)
library(plotly)
p <- ggplot(data = sim_linear, mapping = aes(x = x, y = y)) + 
    geom_point()
ggplotly(p)

1.3 建立回归模型

从散点图上发现:数据点基本排列在一条直线附近，那么我们可以假设\(X\)和\(Y\)的关系是线性的。下面建立以\(zn1.Close\)为自变量，以\(zn2.Close\)为因变量的一元线性模型：可以用公式表式为：

\[Y = \alpha + \beta X + \epsilon\]

• \(Y\) :因变量
• \(X\) :自变量
• \(\alpha\) :截距项
• \(\beta\) :自变量系数
• \(\alpha + \beta X\) :表示\(Y\)随\(X\)的变化而变化的线性部分
• \(\epsilon\) :为残差或随机误差，是其他一切不确定因素影响的总和，其值不可观测。假设\(\epsilon\)是服从均值为0方差为\(\sigma ^ 2\)的正态分布, 记作\(\epsilon \sim N(0, \sigma ^{2})\)

对于上面的公式, 称函数\(f(X)\)为一元线性回归函数:

\[f(X)=\alpha + \beta X\]

• \(\alpha\)为回归常数
• \(\beta\)为回归系数
• \(X\)为回归自变量或回归因子
• \(Y\)为回归因变量或响应变量

如果\((X_1,Y_1), (X_2,Y_2), \ldots, (X_n,Y_n)\)是\((X,Y)\)的一组观测值，则一元线性回归模型可表示为:

\[Y_i = \alpha + \beta X + \epsilon_i, i= 1,2,\ldots,n\]

其中

• \(E(\epsilon_i)=0\)
• \(var(\epsilon_i)=\sigma^2\)
• \(i=1,2,\ldots,n\)

1.4 回归参数估计

简单线性回归模型的参数估计采用最小二乘法，用R语言来实现对上面数据的回归模型的参数估计：

sim_linear_fit <- lm(y ~ 1 + x, data = sim_linear)
# 回归模型的简单信息打印
sim_linear_fit

## 
## Call:
## lm(formula = y ~ 1 + x, data = sim_linear)
## 
## Coefficients:
## (Intercept)            x  
##    118.2874       0.9963

估计出回归参数后，就可以根据回归参数估计得到\(Y\)和\(X\)的一条线性关系直线，称为拟合回归线：

\[Y= \alpha + \beta X = 118.2874 + 0.9963 * X \]

p1 <- p + geom_line(mapping = aes(x = x, y = 118.2874 + 0.9963 * x), colour = "blue")
ggplotly(p1)

p2 <- p + stat_smooth()

拟合回归直线是用数据拟合出来的，是一个近似的值。可以看到有些点在线上，有些点不在线上。要评价这条回归线拟合的好坏，我们就需要对回归模型进行显著性检验。

1.5 回归方程的显著性检验

从回归参数的公式可知，在计算过程中并不一定要知道\(Y\)和\(X\)是否有线性相关的关系。如果不存相关关系，那么回归方程就没有任何意义了，如果\(Y\)和\(X\)是有相关关系的，即\(Y\)会随着\(X\)的变化而线性变化，这个时候一元线性回归方程才有意义。所以，我们需要用假设检验的方法，来验证相关性的有效性。

通常会采用三种显著性检验的方法：

• t检验：t检验是检验模型某个自变量\(X_i\)对于\(Y\)的显著性，即检验：原假设为\(X_i\)的回归系数是否为0。通常用P-value判断显著性，小于0.01或更小时说明这个自变量\(X_i\)与\(Y\)相关关系显著;
• F检验：F检验用于对所有的自变量\(X\)在整体上看对于\(Y\)的线性显著性，即检验：原假设为\(X_i, i= 1,\ldots, p\)的回归系数是否全部为0。也是用P-value判断显著性，小于0.01或更小时说明整体上自变量\(X\)与\(Y\)相关关系显著;
• \(R^2\)(R方)相关系数检验：用来判断回归方程的拟合程度，\(R^2\)的取值在0，1之间，越接近1说明模型对训练数据的拟合程度越好。

在R语言中，上面列出的三种检验的方法都已被实现，我们只需要把结果解读：

# 回归模型的详细信息打印
summary(sim_linear_fit)

## 
## Call:
## lm(formula = y ~ 1 + x, data = sim_linear)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.320 -1.339 -1.283  3.680  8.754 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 118.28740  644.16096   0.184    0.856    
## x             0.99632    0.04563  21.834   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.896 on 25 degrees of freedom
## Multiple R-squared:  0.9502, Adjusted R-squared:  0.9482 
## F-statistic: 476.7 on 1 and 25 DF,  p-value: < 2.2e-16

• Call：列出了回归模型的公式；
• Residuals：列出了残差的最小值点，1/4分位点，中位数点，3/4分位点，最大值点；
• Coefficients: 表示参数估计的计算结果。
  • Estimate: 为参数估计列
    • Intercept: 行表示常数参数的估计值
    • \(x\)行表示自变量x的参数b的估计值。
  • Std. Error: 为参数的标准差，sd(\(\alpha\)), sd(\(\beta\))
  • t value: 为t值，为t检验的值
  • Pr(>|t|): 表示P-value值，用于t检验判定，匹配显著性标记
    • 显著性标记: ***为非常显著, **为高度显著, *为显著，·为不太显著，没有记号为不显著。
• Residual standard error: 表示残差的标准差，自由度为n-2。
• Multiple R-squared: 为相关系数\(R^2\)的检验，越接近1则越显著。
• Adjusted R-squared，为相关系数的修正系数，解决多元回归自变量越多，判定系数R^2越大的问题。
• F-statistic，表示F统计量，自由度为(1,n-2)
  • p-value: 用于F检验判定，匹配显著性标记。

通过查看模型的结果数据，我们可以发现通过T检验的截距和自变量\(X\)都是非常显著，通过F检验判断出整个模型的自变量是非常显著，同时\(R^2\)的相关系数检验可以判断自变量和因变量是高度相关的。

最后，我们通过的回归参数的检验与回归方程的检验，得到最后一元线性回归方程为：

\[Y = 118.2874 + 0.9963 * X\]

1.6 残差分析和异常点检测

在得到的回归模型进行显著性检验后，还要在做残差分析(因变量估计值和实际值之间的差)，检验模型的正确性，残差必须服从正态分布\(\epsilon \sim N(0,\sigma^2)\)。所以需要计算数据残差，并进行正态分布检验。

计算残差：

y.res <- sim_linear$y - fitted(sim_linear_fit)
y.res

##         1         2         3         4         5         6         7 
##  3.550927 -1.375413 -1.375413  3.624587  3.716662 -1.301753 -6.320168 
##         8         9        10        11        12        13        14 
## -1.320168 -1.338583 -1.338583 -6.283338 -1.338583 -1.338583 -1.320168 
##        15        16        17        18        19        20        21 
##  3.661417 -1.283338 -1.264923  3.698247 -6.228092 -1.264923  3.771908 
##        22        23        24        25        26        27 
##  3.790323 -6.209677  3.771908 -1.209677  3.771908  8.753493

y.res <- residuals(sim_linear_fit)
y.res

##         1         2         3         4         5         6         7 
##  3.550927 -1.375413 -1.375413  3.624587  3.716662 -1.301753 -6.320168 
##         8         9        10        11        12        13        14 
## -1.320168 -1.338583 -1.338583 -6.283338 -1.338583 -1.338583 -1.320168 
##        15        16        17        18        19        20        21 
##  3.661417 -1.283338 -1.264923  3.698247 -6.228092 -1.264923  3.771908 
##        22        23        24        25        26        27 
##  3.790323 -6.209677  3.771908 -1.209677  3.771908  8.753493

参差正态性检验：

shapiro.test(y.res)

## 
##  Shapiro-Wilk normality test
## 
## data:  y.res
## W = 0.86459, p-value = 0.002272

res = ggplot(data = data.frame(cbind(index = 1:length(y.res)), y.res),
             mapping = aes(x = index, y = y.res)) +
    geom_point()
ggplotly(res)

对残差进行Shapiro-Wilk正态分布检验，W接近1，p-value>0.05，证明数据集符合正态分布！

生成评价拟合模型拟合情况的四幅诊断图：

par(mfrow = c(2, 2))
plot(sim_linear_fit)

• 图1:残差和拟合值对比图
  • 对残差和拟合值作图，横坐标是拟合值，纵坐标是残差。残差和拟合值之间，数据点均匀分布在y=0两侧，呈现出随机的分布，红色线呈现出一条平稳的曲线并没有明显的形状特征，说明残差数据表现非常好。
• 图2:残差QQ图
  • 残差QQ图，用来描述残差是否符合正态分布。图中的数据点按对角直线排列，趋于一条直线，并被对角直接穿过，直观上符合正态分布。对于近似服从正态分布的标准化残差，应该有 95% 的样本点落在 [-2,2] 区间内。
• 图3:标准化残差平方根和拟合值对比图
  • 对标准化残差平方根和拟合值作图，横坐标是拟合值，纵坐标是标准化后的残差平方根。与残差和拟合值对比图(图1)的判断方法类似，数据随机分布，红色线呈现出一条平稳的曲线，无明显的形状特征。
• 图4:标准残差和杠杆值对比图
  • 对标准化残差和杠杆值作图，虚线表示的cooks距离等高线，通常用Cook距离度量的回归影响点。本图中没有出现红色的等高线，则说明数据中没有特别影响回归结果的异常点。

看到上面4幅中，每幅图上都有一些点被特别的标记出来了，这些点是可能存在的异常值点，如果要对模型进行优化，我们可以从这些来入手。但终于本次残差分析的结果已经很好了，所以对于异常点的优化，可能并不能明显的提升模型的效果。

1.7 模型预测

通过上面的建模，获得了一元线性回归方程的公式，就可以对数据进行预测了。

对给定\(X=x_0\)时，计算出\(y_0=\alpha +\beta*x_0\)的值，并计算出置信度为\(1-\alpha\)的预测区间。

当\(X=x_0\),\(Y=y_0\)时，置信度为\(1-\alpha\)的预测区间为:

\[[\hat{y}_0-l, \hat{y}_0+l]\]

其中：

\[l=t_{\alpha}(n-2)\hat{\sigma}\sqrt{1+\frac{1}{n}+\frac{(\bar{x}-x_0)^2}{S_{xx}}}\]

即：

\[P(\hat{y}_0-l <y_{0} <\hat{y}_0+l) = \alpha\]

计算预测值y0，和相应的预测区间

new_x <- data.frame(x=14113)
sim_linear_fit_pred <- predict(sim_linear_fit, new_x, interval = "prediction", level = 0.95)
sim_linear_fit_pred

##        fit      lwr      upr
## 1 14179.31 14171.13 14187.49

当\(x_0=14113\)时，在预测区间为0.95的概率时，\(y_0\)的值为14106.58，预测区间为\([14095.69, 14117.46]\)。

我们通过图形来表示。

p <- ggplot(data = sim_linear, mapping = aes(x = x, y = y)) + 
    geom_point() + 
    geom_line(mapping = aes(x = x, y = 118.2874 + 0.9963 * x),  colour = "yellow") + 
    geom_point(data = data.frame(x1 = new_x[, 1], 
                                 y1 = sim_linear_fit_pred[1:3]),
               mapping = aes(x = x1, y = y1), 
               colour = c("red", "blue", "green"))
ggplotly(p)

其中，红色点为\(y_0\)的值，蓝色点为预测区间最小值，绿色点为预测区间最大值。