77-预测分析－Ｒ语言实现-线性回归
2020-09-25 本文已影响0人
wonphen
１、线性回归假设
假设：
１、假设输出变量是一组特征变量的加权线性函数；
２、对于特征变量的固定值，输出是具有常数方差的正态分布；
３、特征本身在统计学上是相互独立的。
协方差衡量的是两个变量之间关联的紧密程度，可以是正值或负值。正的协方差表明了正相关性，也就是说，当一个变量增大时，另一个也会增大。负的协方差则代表相反的含义，当一个变量增大时，另一个往往会减小。当两个变量从统计学角度是互相独立因而不相关时，它们的协方差就是0，反之则不一定。
２、导入数据
library(pacman)
p_load(dplyr, readr, caret)
machine <- read_csv("./data_set/machine.data", col_names = F)
str(machine)
## tibble [209 × 10] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ X1 : chr [1:209] "adviser" "amdahl" "amdahl" "amdahl" ...
## $ X2 : chr [1:209] "32/60" "470v/7" "470v/7a" "470v/7b" ...
## $ X3 : num [1:209] 125 29 29 29 29 26 23 23 23 23 ...
## $ X4 : num [1:209] 256 8000 8000 8000 8000 8000 16000 16000 16000 32000 ...
## $ X5 : num [1:209] 6000 32000 32000 32000 16000 32000 32000 32000 64000 64000 ...
## $ X6 : num [1:209] 256 32 32 32 32 64 64 64 64 128 ...
## $ X7 : num [1:209] 16 8 8 8 8 8 16 16 16 32 ...
## $ X8 : num [1:209] 128 32 32 32 16 32 32 32 32 64 ...
## $ X9 : num [1:209] 198 269 220 172 132 ...
## $ X10: num [1:209] 199 253 253 253 132 ...
## - attr(*, "spec")=
## .. cols(
## .. X1 = col_character(),
## .. X2 = col_character(),
## .. X3 = col_double(),
## .. X4 = col_double(),
## .. X5 = col_double(),
## .. X6 = col_double(),
## .. X7 = col_double(),
## .. X8 = col_double(),
## .. X9 = col_double(),
## .. X10 = col_double()
## .. )
# 添加列名
names(machine) <- c("vendor", "model", "myct", "mmin", "mmax",
"cach", "chmin", "chmax", "prp", "erp")
machine <- machine[, 3:9]
dim(machine)
## [1] 209 7
DataExplorer::profile_missing(machine)
## # A tibble: 7 x 3
## feature num_missing pct_missing
## <fct> <int> <dbl>
## 1 myct 0 0
## 2 mmin 0 0
## 3 mmax 0 0
## 4 cach 0 0
## 5 chmin 0 0
## 6 chmax 0 0
## 7 prp 0 0
数据集没有缺失值。
３、拆分训练集和测试集
set.seed(123)
ind <- createDataPartition(machine$prp, p = 0.85, list = F)
dtrain <- machine[ind, ]
dtest <- machine[-ind, ]
dim(dtrain)
## [1] 179 7
dim(dtest)
## [1] 30 7
４、设置预处理参数
设置相关性截断阈值为0.75。
ct <- trainControl(preProcOptions = list(cutoff = 0.75))
５、建立线性回归模型
预处理函数为剔除高相关性函数。
fit <- train(prp ~ ., data = dtrain, method = "lm",
trControl = ct, preProcess = c("corr"))
summary(fit$finalModel)
##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -168.64 -28.93 4.72 28.45 345.00
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.916e+01 9.002e+00 -6.572 5.70e-10 ***
## myct 5.425e-02 2.005e-02 2.706 0.0075 **
## mmin 1.528e-02 2.095e-03 7.294 1.06e-11 ***
## mmax 5.722e-03 6.924e-04 8.265 3.63e-14 ***
## cach 5.992e-01 1.485e-01 4.035 8.20e-05 ***
## chmin -9.309e-01 1.074e+00 -0.867 0.3871
## chmax 1.753e+00 2.442e-01 7.180 2.00e-11 ***
## ---
## Signif. codes:
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 61.37 on 172 degrees of freedom
## Multiple R-squared: 0.865, Adjusted R-squared: 0.8603
## F-statistic: 183.7 on 6 and 172 DF, p-value: < 2.2e-16
６、残差分析
# 残差分布
summary(fit$finalModel$residuals)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -168.643 -28.932 4.719 0.000 28.449 344.997
# 输出变量的平均值
mean(dtrain$prp)
## [1] 104.9553
对比残差的值和输出变量的平均值，以查看残差是否过大。预测结果的50%完全都在正确值+-2９范围内，结果相当合理。
创建残差分位图：
qqnorm(fit$finalModel$residuals, main = "Ｎormal Q-Q Plot for CPU data set")
qqline(fit$finalModel$residuals)
残差分位图
残差看上去跟理论上的正态分布分位数相当贴近，拟合程度并非完美，但是这对于大部分实际环境来说是典型的情况。
# 另一种获得残差图的方法
par(mfrow = c(2, 2))
plot(fit$finalModel)
残差图
查看拟合值－残差图（第一张）可知，左侧部分显示了细微的残差递减模式，这表明残差可能并不是同方差的（违背线性回归假设）。
７、计算Ｐ值
手动计算chmin回归系数的Ｐ值，因为为双边检验，所以乘以２，否则的话就是单边检验。
自由度df的数量是在计算某个统计量（例如某个系数的估计值）时能够自由变化的变量数。在线性回归中，它就相当于训练数据中观测数据的个数减去模型中参数的个数。自由度这个名字来源于它和独立维度的数量之间的关系，从而反应系统能够在不违反任何对于输入的约束条件的条件下自由配置的程度。
对Ｐ值的理解：
１、不能通过比较Ｐ值来评判哪个特征更重要；
２、较高的Ｐ值并不一定代表某个特征和输出之间不存在线性关系，它只表明在所有其他特征存在时，这个特征对于输出不能提供任何新的信息；
３、９５％的经验并非绝对真理，当特征数不是特别多时有用，但对于高维问题并不是那么有用。
t.value <- -9.309e-01 / 1.074e+00
p.value <- pt(t.value, df = 172) * 2
p.value
## [1] 0.3872812
与summary中计算的结果一致。
８、方差分析
F统计量是用来衡量两个正态分布的方差之间是否存在统计显著性，本实例中的Ｆ统计量用来评价所有系数都为０的模型中的残差方差是否和我们训练的模型的残差方差有显著性差异。
fit.null <- lm(prp ~ 1, data = dtrain)
fit.full <- lm(prp ~ ., data = dtrain)
anova(fit.null, fit.full)
## Analysis of Variance Table
##
## Model 1: prp ~ 1
## Model 2: prp ~ myct + mmin + mmax + cach + chmin + chmax
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 178 4798642
## 2 172 647806 6 4150835 183.68 < 2.2e-16 ***
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
与summary中计算的一致：F-statistic: 183.7 on 6 and 172 DF, p-value: < 2.2e-16
９、检测多重共线性
多重共线性表现：
１ 两个高度共线性的特征具有较大的Ｐ值，表明它们和输出变量不相关，但如果去除其中一个特征并重新训练模型，剩下的那个特征就会有较小的Ｐ值；
２ 某个系数出现不正常的符号。
car::vif(fit$finalModel)
## myct mmin mmax cach chmin chmax
## 1.209972 3.070319 3.254166 1.762142 2.039244 1.903203
VIF分数为４或更大的特征是可疑的，分数大于１０就表明多重共线性的极大可能性．本实例全部小于４，说明不存在多重共线性问题。
１０、模型评价
均方误差ＭＳＥ。
compute_mse <- function(prediction, actual) {
mean((prediction - actual) ^ 2)
}
compute_mse(fit.full$fitted.values, dtrain$prp)
## [1] 3619.029
11、离群值
从残差分析图可以看到Ｘ１７３行观测出现在了四张图中，说明它极有可能是离群值．去掉该观测后重新训练模型试试。
p_load(broom)
fit.omit <- lm(prp ~ ., data = dtrain[!(rownames(dtrain)) %in% c(173), ])
compute_mse(fit.omit$fitted.values, dtrain[!(rownames(dtrain)) %in% c(173), ]$prp)
## [1] 2690.545
data.frame(fit.full = t(glance(fit$finalModel))),
fit.omit = t(glance(fit.omit)))
## fit.full fit.omit
## r.squared 8.650022e-01 8.705765e-01
## adj.r.squared 8.602930e-01 8.660354e-01
## sigma 6.137031e+01 5.292149e+01
## statistic 1.836824e+02 1.917073e+02
## p.value 4.641565e-72 3.437299e-73
## df 6.000000e+00 6.000000e+00
## logLik -9.873495e+02 -9.554485e+02
## AIC 1.990699e+03 1.926897e+03
## BIC 2.016198e+03 1.952351e+03
## deviance 6.478061e+05 4.789170e+05
## df.residual 1.720000e+02 1.710000e+02
## nobs 1.790000e+02 1.780000e+02
去掉离群值后R2略微增加了，ＭＳＥ也降低了不少。
12、预测
caret.pred <- predict(fit, newdata = dtest)
omit.pred <- predict(fit.omit, newdata = dtest)
data.frame(caret_mse = compute_mse(caret.pred, dtest$prp),
omit_mse = compute_mse(omit.pred, dtest$prp))
## caret_mse omit_mse
## <dbl> <dbl>
## 2914.014 1672.859
去掉离群值后，模型在测试集上的ＭＳＥ也小了很多．
