多元线性回归实战模型建立

个人总结的部分多元线性回归案例，目前线性回归过于简单，在机器学习盛行的当下在数学建模和数据挖掘里面应用较少，但是该模型简单容易上手。

多元线性回归模型表达式：

多元线性回归，顾名思义肯定是解决多个自变量变量对一个因变量的线性问题

案例1：美国房价预测

下面以美国房价案例来进行分析

##多元线性回归
####案例1：美国房屋信息数据的预测####
##预测AvgPrice，并建立多元线性回归模型
##USA_Housing.csv包含不同地区的平均房价及多个可能影响房价的自变量（PS：我个人翻译可能会有问题）
##AvgAreaIncome：该地区的平均收入
##AvgAreaHouseAge：房子的平均面积
##AvgAreaNumberRooms：房子平均数量
##AvgAreaNumberofBedrooms：平均卧室面积
##AreaPopulation：该地区人口数量
##AvgPrice：房子价格
getwd()##寻找默认目录，将USA_Housing.csv数据集放到该目录下
data <- read.csv("USA_Housing.csv")##读取数据集
summary(data)##观察数据，存在6个变量

##作密度图，可视化各个变量的分布
library(ggcorrplot)
library(tidyr)
data_0 <- gather(data,key="varname",value="value",1:6)
ggplot(data_0)+theme_bw()+
  geom_density(aes(value),fill = "red",alpha = 0.5)+
  facet_wrap(.~varname,scales = "free")+
  theme(axis.text.x = element_text(angle = 30))
##除了AvgAreaNumberofBedrooms以外其它变量都呈现正态分布

除了AvgAreaNumberofBedrooms以外其它变量都呈现正态分布，一般来说，多元线性回归的前提条件需要保证自变量近似正态分布,。

##作热力图观察变量之间的关联性
house_cor <- cor(data)
library(corrplot)##作热力图必须的包
corrplot.mixed(house_cor,tl.col="black",tl.pos = "d",number.cex = 0.8)

##换一种热力图作法
library(ggcorrplot)
ggcorrplot(house_cor,method = "square",lab = TRUE)+
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10))
##颜色越深，表面变量之间的关联性越强，AvgPrice除了和AvgAreaNumberofBedrooms，
##都明显呈现出关联性，且和AvgAreaIncome关联性最强=0.64,说明可以建立多元线性回归模型。

颜色越深，表面变量之间的关联性越强，AvgPrice除了和AvgAreaNumberofBedrooms，都明显呈现出关联性，且和AvgAreaIncome关联性最强=0.64,说明可以建立多元线性回归模型。

##多元线型回归的建立
lm1 <- lm(AvgPrice~.,data = data)
summary(lm1)
##一般认为，P值小于0.05则认为变量显著，从输出的结果中可以发现AvgAreaNumberofBedrooms变量是不显著的p值等于0.207
##但是整体模型的p值小于0.05，且R2等于0.9179，建立得到的模型是显著的

第一次回归结果：

Call:
lm(formula = AvgPrice ~ ., data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-337020  -70236     320   69175  361870 

Coefficients:
                          Estimate Std. Error  t value Pr(>|t|)    
(Intercept)             -2.637e+06  1.716e+04 -153.708   <2e-16 ***
AvgAreaIncome            2.158e+01  1.343e-01  160.656   <2e-16 ***
AvgAreaHouseAge          1.656e+05  1.443e+03  114.754   <2e-16 ***
AvgAreaNumberRooms       1.207e+05  1.605e+03   75.170   <2e-16 ***
AvgAreaNumberofBedrooms  1.651e+03  1.309e+03    1.262    0.207    
AreaPopulation           1.520e+01  1.442e-01  105.393   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 101200 on 4994 degrees of freedom
Multiple R-squared:  0.918,	Adjusted R-squared:  0.9179 
F-statistic: 1.119e+04 on 5 and 4994 DF,  p-value: < 2.2e-16

##这里，为了提高模型精度，剔除不显著变量拟合新的回归模型
lm2 <- lm(AvgPrice~AvgAreaIncome+AvgAreaHouseAge+AvgAreaNumberRooms
          +AreaPopulation,data = data)
summary(lm2)
##剔除AvgAreaNumberofBedrooms后，模型p值没有变化，R2=0.9179，说明剔除该变量对模型没有影响。

剔除变量第二次回归结果：

Call:
lm(formula = AvgPrice ~ AvgAreaIncome + AvgAreaHouseAge + AvgAreaNumberRooms + 
    AreaPopulation, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-338419  -70058     132   69074  362025 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -2.638e+06  1.716e+04 -153.73   <2e-16 ***
AvgAreaIncome       2.158e+01  1.343e-01  160.74   <2e-16 ***
AvgAreaHouseAge     1.657e+05  1.443e+03  114.77   <2e-16 ***
AvgAreaNumberRooms  1.216e+05  1.423e+03   85.48   <2e-16 ***
AreaPopulation      1.520e+01  1.442e-01  105.39   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 101200 on 4995 degrees of freedom
Multiple R-squared:  0.918,	Adjusted R-squared:  0.9179 
F-statistic: 1.398e+04 on 4 and 4995 DF,  p-value: < 2.2e-16

后面为多元线性回归检验步骤，一般在数学建模里面，不单单要建立模型，还得对模型进行检验与评估，分析模型的好坏。

##多元线性回归模型的检验
##可视化回归模型的系数
library(GGally)
ggcoef(lm2,exclude_intercept = T,vline_color = "red",
       errorbar_color = "blue",errorbar_height = 0.1)+
  theme_bw()
## AvgAreaNumberRooms和AvgAreaHouseAge的回归系数最大，说明他们和AvgPrice关联性高

## 可视化回归模型的图像
par(mfrow = c(2,2))
plot(lm2)
## 图中四个图从左到右分别表示：残差-拟合图(Residuals vs Fitted)、标准化残差的正态性检验Q-Q图(NormalQ-Q)、
## 标准化残差-拟合图(Scale-Local)、残差一杠杆图(Residuals vs Leverage)。从4个子图中可以发现，除几个特殊样
## 本(4717、40等)，模型拟合的残差属于正态分布，即该数据集满足多元线性回归的条件，且回归的结果是显著的。

## 对此可以建立多元线性回归的表达式
## AvgPrice=-2638000+21.58*AvgArealncome+165700*AvgAreaHouseAge+121600*AvgAreaNumberRooms+15.2*AreaPopulation

最后得到多元线性回归的表达式：
AvgPrice=-2638000+21.58*AvgArealncome+165700*AvgAreaHouseAge+121600*AvgAreaNumberRooms+15.2*AreaPopulation
 

案例2：囊胞纤维症患者肺功能的研究

下列为囊胞纤维症患者肺功能数据

以fev1为因变量，建立预测模型。

##cystfibr是ISwR包中的数据集，该数据帧具有25行和10列。它包含囊性纤维化患者（7–23岁）的肺功能数据。
##age，年龄
##sex，性别
##height，身高，单位:cm
##weight，体重，单位:kg
##bmp，与正常体重的百分体重，单位:%
##fev1，最大肺活量
##rv，余气量
##frc，功能残气量
##tlc，肺总气量
##pemax，最大呼气压力
##（PS:看不懂的上面名词话建议学一学人体生理解剖学）
##建立fev1（最大肺活量）的预测模型。
library(ISwR)
cystfibr<-cystfibr
summary(cystfibr)

##作密度图，可视化各个变量的分布
library(ggcorrplot)
library(tidyr)
data_1 <- gather(cystfibr,key="varname",value="value",1:10)
ggplot(data_1)+theme_bw()+
  geom_density(aes(value),fill = "red",alpha = 0.5)+
  facet_wrap(.~varname,scales = "free")+
  theme(axis.text.x = element_text(angle = 30))
##除bmp和sex以外，其它变量都近似呈正态分布。

除bmp和sex以外，其它变量都近似呈正态分布。

##作热力图观察变量之间的关联性
cystfibr_cor <- cor(cystfibr)
library(ggcorrplot)
ggcorrplot(cystfibr_cor,method = "square",lab = TRUE)+
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10))
##不难观察得到，自变量和因变量fev1存在明显的相关性，大部分绝对值都大于0.4以上。且frc关联性最强，为-0.67。
##但是，观察得到自变量rv,frc,tlc之间存在明显的关联性，说明可能存在多重共线性问题。

观察得到自变量rv,frc,tlc之间存在明显的关联性，说明可能存在多重共线性问题。

##先不讨论多重共线性问题，以因变量fev1直接建立多元线性回归模型
fy1<-lm(fev1~.,data = cystfibr)
summary(fy1)
##最后得到模型的R2等于0.7232，p值等于0.0002646<0.05，说明多元线性效果还是可以的。
##但是，观察自变量的p值不难发现，只有4个自变量的p值小于0.05是显著的，说明其它自变量对提高模型稳定性可能没有帮助。

Call:
lm(formula = fev1 ~ ., data = cystfibr)

Residuals:
   Min     1Q Median     3Q    Max 
-7.278 -4.113  1.565  3.052  8.126 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept) 104.98313   45.88316   2.288  0.03707 * 
age          -1.90901    1.00668  -1.896  0.07734 . 
sex          -9.24703    2.67096  -3.462  0.00348 **
height       -0.24728    0.20071  -1.232  0.23690   
weight        0.65133    0.46833   1.391  0.18459   
bmp           0.03319    0.28666   0.116  0.90937   
rv            0.04594    0.04536   1.013  0.32722   
frc          -0.24203    0.09697  -2.496  0.02470 * 
tlc          -0.07717    0.11440  -0.675  0.51022   
pemax         0.05781    0.05782   1.000  0.33328   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.891 on 15 degrees of freedom
Multiple R-squared:  0.827,	Adjusted R-squared:  0.7232 
F-statistic: 7.968 on 9 and 15 DF,  p-value: 0.0002646

多元线性回归模型的检验

##可视化回归模型的系数
library(GGally)
ggcoef(fy1,exclude_intercept = T,vline_color = "red",
       errorbar_color = "blue",errorbar_height = 0.1)+
  theme_bw()
## AvgAreaNumberRooms和AvgAreaHouseAge的回归系数最大，说明他们和AvgPrice关联性高

## 可视化回归模型的图像
par(mfrow = c(2,2))
plot(fy1)
## Q-Q图不难观察得到，除了个别样本（4，10，18），模型拟合的残差属于正态分布。

多重共线性问题与逐步回归

回归通常都会出现两种问题，过拟合和欠拟合，其中多重共线性的问题与过拟合密切相关。过拟合通常会使得模型失去泛用性，欠拟合会导致模型拟合的精度下降，前者可以通过减少无关变量来提高模型效果。继续分析案例2
案例2中最后的模型发现存在多重共线性问题

cystfibr_cor##查看相关性系数
##不难观察得到自变量很多之间都存在明显的相关性，age、height、weight之间相关性系数都大于0.9

##计算模型的方程膨胀因子（VIF），VIF值通常表示自变量之间共线性程度。
data_VIF<-cystfibr[,c(1:5,7:10)]##注意，计算自变量之间的VIF值要剔除掉因变量。
fy1_VIF<-summary(lm(age~.,data = data_VIF))##这里设置年龄为因变量
1/(1-fy1_VIF$r.squared)
##计算得到VIF值等于17.9，一般情况认为VIF值大于5时，则可以认为存在多重共线性问题。

##vif()是car包中特有的可以计算模型全部自变量的VIF函数。
library(car)
vif(fy1)##观察全部自变量的VIF值

逐步回归模型的建立

##逐步回归##
## 逐步回归通常用于优化多变量的线性回归模型，它能够剔除不显著变量建立回归模型。
## 赤池信息准则（AIC），该准则用于分析模型的拟合优度和用来拟合参数的数目，
## AIC值越小说明多重共线性小。
drop1(fy1)##观察模型的AIC值
##结果观察得到，当不剔除任何变量时，AIC值等于95。
drop1(fy1,test = "F")
## 分析得到最适合剔除的变量为sex
## 在R语言中，运用step()函数自动选择最优的方法进行多元线性回归模型。对于R语言
## 参数direction可以设置为"both"(向前向后法)、"backward"(向后法)或"forward"(向前法)。
fy2<-step(fy1)
## 结果表明只使用height、rv、age、weight、frc和sex为自变量进行回归的AIC值最小。
summary(fy2)
## 逐步回归前模型：R2 = 0.7232，p值 = 0.0002646
## 逐步回归后模型：R2 = 0.7481，p值 = 1.163e-05
## 说明使用逐步回归对模型拟合的效果是显著的，同时模型的拟合度和置信度都有所提高。

Call:
lm(formula = fev1 ~ age + sex + height + weight + rv + frc, data = cystfibr)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.3378 -3.8380  0.4663  3.9743  7.3457 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 108.34533   19.46710   5.566 2.78e-05 ***
age          -2.02046    0.70879  -2.851 0.010619 *  
sex          -9.96707    2.39768  -4.157 0.000592 ***
height       -0.26152    0.16556  -1.580 0.131600    
weight        0.77778    0.20089   3.872 0.001118 ** 
rv            0.06688    0.03876   1.725 0.101572    
frc          -0.30219    0.07709  -3.920 0.001004 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.62 on 18 degrees of freedom
Multiple R-squared:  0.8111,	Adjusted R-squared:  0.7481 
F-statistic: 12.88 on 6 and 18 DF,  p-value: 1.163e-05

逐步回归模型的检验

##模型的检验
confint(fy2)## t检验
anova(fy2)## F检验
##结果表面rv、age、weight、frc和sex为自变量进行回归结果是显著的。

##回归诊断
library(gvlma)
gvlma(fy2)
##结果全局统计量(Global Stat)与其他个检验项的p值均大于0.05，数据满足线性回归所以的统计假设。
##epiDisplay包中的regress.display()函数可以汇总模型的结果。
library(epiDisplay)
regress.display(fy2)

模型检验结果：
Linear regression predicting fev1
 
                    adj. coeff.(95%CI)   P(t-test) P(F-test)
age (cont. var.)    -2.02 (-3.51,-0.53)  0.011     0.01     
                                                            
sex: 1 vs 0         -9.97 (-15,-4.93)    < 0.001   < 0.001  
                                                            
height (cont. var.) -0.26 (-0.61,0.09)   0.132     0.343    
                                                            
weight (cont. var.) 0.78 (0.36,1.2)      0.001     0.001    
                                                            
rv (cont. var.)     0.07 (-0.01,0.15)    0.102     < 0.001  
                                                            
frc (cont. var.)    -0.3 (-0.46,-0.14)   0.001     0.001    
                                                            
No. of observations = 25