R语言学习系列32-回归分析(共26页).docx

精选优质文档-倾情为你奉上27. 回归分析回归分析是研究一个或多个变量（因变量）与另一些变量（自变量）之间关系的统计方法。主要思想是用最小二乘法原理拟合因变量与自变量间的最佳回归模型（得到确定的表达式关系）。其作用是对因变量做解释、控制、或预测。回归与拟合的区别：拟合侧重于调整曲线的参数，使得与数据相符；而回归重在研究两个变量或多个变量之间的关系。它可以用拟合的手法来研究两个变量的关系，以及出现的误差。 回归分析的步骤： （1）获取自变量和因变量的观测值；（2）绘制散点图，并对异常数据做修正；（3）写出带未知参数的回归方程；（4）确定回归方程中参数值；（5）假设检验，判断回归方程的拟合优度；（6）进行解释、控制、或预测。（一）一元线性回归一、原理概述 1. 一元线性回归模型： Y=𝛽0+𝛽1X+其中 X是自变量，Y是因变量，𝛽0,𝛽1是待求的未知参数，𝛽0也称为截距；是随机误差项，也称为残差，通常要求满足： 的均值为0； 的方差为𝜎2； 协方差COV(i, j)=0，当ij时。即对所有的ij, i与j 互不相关。用最小二乘法原理，得到最佳拟合效果的值：， 2.模型检验(1) 拟合优度检验计算R2，反映了自变量所能解释的方差占总方差的百分比，值越大说明模型拟合效果越好。通常可以认为当R2大于0.9时，所得到的回归直线拟合得较好，而当R2小于0.5时，所得到的回归直线很难说明变量之间的依赖关系。(2) 回归方程参数的检验回归方程反应了因变量Y随自变量X变化而变化的规律，若𝛽1=0，则Y不随X变化，此时回归方程无意义。所以，要做如下假设检验：H0: 𝛽1=0, H1: 𝛽10； F检验 若𝛽1=0为真，则回归平方和RSS与残差平方和ESS/(N-2)都是𝜎2的无偏估计，因而采用F统计量：来检验原假设1=0是否为真。 T检验对H0: 𝛽1=0的T检验与F检验是等价的（t2=F）。3. 用回归方程做预测得到回归方程后，预测X=x0处的Y值.的预测区间为：其中t/2的自由度为N-2. 二、R语言实现使用lm()函数实现，基本格式为：lm(formula, data, subset, weights, na.action,method="qr", .)其中，formula为要拟合的回归模型的形式，一元线性回归的格式为：yx，y表示因变量，x表示自变量，若不想包含截距项，使用yx-1；data为数据框或列表；subset选取部分子集；weights取NULL时表示最小二乘法拟合，若取值为权重向量，则用加权最小二乘法；na.action设定是否忽略缺失值；method指定拟合的方法，目前只支持“qr”（QR分解），method=“model.frame”返回模型框架。三、实例例1 现有埃及卡拉马村庄每月记录儿童身高的数据，做一元线性回归。datas<-data.frame(age=18:29,height=c(76.1,77,78.1,78.2,78.8,79.7,79.9,81.1,81.2,81.8,82.8,83.5)datas age height1 18 76.12 19 77.03 20 78.14 21 78.25 22 78.86 23 79.77 24 79.98 25 81.19 26 81.210 27 81.811 28 82.812 29 83.5plot(datas) #绘制散点图res.reg<-lm(heightage,datas) #做一元线性回归summary(res.reg) #输出模型的汇总结果Residuals: Min 1Q Median 3Q Max -0.27238 -0.24248 -0.02762 0.16014 0.47238 Coefficients: Estimate Std.Error t value Pr(>|t|) (Intercept) 64.9283 0.5084 127.71 < 2e-16 *age 0.6350 0.0214 29.66 4.43e-11 *-Signif.codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Residual standard error: 0.256 on 10 degrees of freedomMultiple R-squared: 0.9888,Adjusted R-squared: 0.9876 F-statistic: 880 on 1 and 10 DF, p-value: 4.428e-11说明：输出了残差信息Residuals；回归系数估计值、标准误、t统计量值、p值，可得到回归方程：height=64.9283+0.6350*age回归系数p值（< 2e-16，4.43e-11）很小，非常显著的0；*也表示显著程度非常显著。拟合优度R2=0.9888>0.5, 表示拟合程度很好。F统计量=880，p值=4.428e-11远小于0.05，表示整个回归模型显著，适合估计height这一因变量。coefficients(res.reg) #返回模型的回归系数估计值(Intercept) age 64. 0.confint(res.reg,parm="age",level=0.95) #输出参数age的置信区间，若不指定parm将返回所有参数的置信区间 2.5 % 97.5 %age 0. 0.fitted(res.reg) #输出回归模型的预测值 1 2 3 4 5 6 7 8 9 10 11 1276.35769 76.99266 77.62762 78.26259 78.89755 79.53252 80.16748 80.80245 81.43741 82.07238 82.70734 83.34231anova(res.reg) #输出模型的方差分析表Response: height Df Sum Sq Mean Sq F value Pr(>F) age 1 57.655 57.655 879.99 4.428e-11 *Residuals 10 0.655 0.066 -Signif.codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1vcov(res.reg) #输出模型的协方差矩阵 (Intercept) age(Intercept) 0. -0.age -0. 0.residuals(res.reg) #输出模型的残差 1 2 3 4 5 6 7 8 9 10 11 12-0. 0. 0. -0. -0. 0. -0. 0. -0. -0. 0. 0.AIC(res.reg) #输出模型的AIC值1 5.BIC(res.reg) #输出模型的BIC值1 6.logLik(res.reg) #输出模型的对数似然值'log Lik.' 0. (df=3)abline(res.reg) #给散点图加上一条回归线par(mfrow=c(2,2)plot(res.reg) #绘制回归诊断图说明：分别是残差与拟合值图，二者越无关联越好，若有明显的曲线关系，则说明需要对线性回归模型加上高次项；残差的Q-Q图，看是否服从正态分布；标准化残差与拟合值图，也叫位置-尺度图，纵坐标是标准化残差的平方根，残差越大，点的位置越高，用来判断模型残差是否等方差，若满足则水平线周围的点应随机分布；残差与杠杆图，虚线表示Cooks距离（每个数据点对回归线的影响力）等高线，从中可以鉴别出离群点（第3个点较大，表示删除该数据点，回归系数将有实质上的改变，为异常值点）、高杠杆点、强影响点。datas<-datas-3, #删除第3个样本点，重新做一元线性回归res.reg2<-lm(heightage,datas)summary(res.reg2) 新的回归方程为：height=64.5540+0.6489*age，拟合优度R2=0.993，拟合效果变得更好。 #用回归模型预测ages<-data.frame(age=30:34)pre.res<-predict(res.reg2,ages,interval="prediction",level=0.95) #注意predict函数的第1个参数必须是回归模型的自变量数据构成的数据框或列表pre.res fit lwr upr1 84.02034 83.46839 84.572282 84.66921 84.09711 85.241323 85.31809 84.72365 85.912544 85.96697 85.34825 86.585695 86.61585 85.97114 87.26056（二）多元线性回归一、基本原理 1. 多元线性回归模型：Y=𝛽0+𝛽1X1+ 𝛽NXN+其中 X1, , XN是自变量，Y是因变量，𝛽0, 𝛽1, 𝛽N是待求的未知参数，是随机误差项（残差），若记多元线性回归模型可写为矩阵形式：Y=X+通常要求：矩阵X的秩为k+1（保证不出现共线性）, 且k<N; 为正态分布，E()=0 和E()= 𝜎2I，其中I为N×N单位矩阵。用最小二乘法原理，令残差平方和 最小，得到为的最佳线性无偏估计量（高斯马尔可夫定理）。2. 𝜎2的估计和T检验选取𝜎2的估计量：则假如t值的绝对值相当大，就可以在适当选定的置信水平上否定原假设，参数的1-置信区间可由下式得出：其中t/2为与%显著水平有关的t分布临界值。3. R2和F检验若因变量不具有0平均值，则必须对R2做如下改进：随着模型中增添新的变量，R2的值必定会增大，为了去掉这种增大的干扰，还需要对R2进行修正（校正拟合优度对自由度的依赖关系）：做假设检验：H0: 𝛽1=𝛽N=0; H1: 𝛽1, 𝛽N至少有一个0；使用F统计量做检验，若F值较大，则否定原假设。4. 回归诊断（1）残差图分析残差图就是以残差为纵坐标，某一个合适的自变量为横坐标的散点图。回归模型中总是假定误差项是独立的正态分布随机变量，且均值为零和方差相等为𝜎2. 如果模型适合于观察到的数据，那么残差作为误差的无偏估计，应基本反映误差的假设特征。即残差图应该在零点附近对称地密布，越远离零点的地方就疏散（在形象上似有正态趋势），则认为模型与数据拟合得很好。若残差图呈现如图（a）所示的形式，则认为建立的回归模型正确，更进一步再诊断“学生化残差”是否具有正态性：图（b）表明数据有异常点，应处理掉它重新做回归分析（在SAS的REG回归过程步中用来度量异常点影响大小的统计量是COOKD统计量）；图（c）残差随x的增大而增大，图（d）残差随x的增大而先增后减，都属于异方差。此时应该考虑在回归之前对数据y或x进行变换，实现方差稳定后再拟合回归模型。原则上，当误差方差变化不太快时取变换；当误差方差变化较快时取变换log y或ln y；当误差方差变化很快时取变换1/y；还有其他变换，如著名的Box-Cox幂变换.图（e）（f）表示选用回归模型是错误的。（2）共线性回归分析中很容易发生模型中两个或两个以上的自变量高度相关，从而引起最小二乘估计可能很不精确（称为共线性问题）。在实际中最常见的问题是一些重要的自变量很可能由于在假设检验中t值不显著而被不恰当地剔除了。共线性诊断问题就是要找出哪些变量间存在共线性关系。（3）误差的独立性回归分析之前，要检验误差的独立性。若误差项不独立，那么回归模型的许多处理，包括误差项估计、假设检验等都将没有推导依据。由于残差是误差的合理估计，因此检验统计量通常是建立在残差的基础上。检验误差独立性的最常用方法，是对残差的一阶自相关性进行Durbin-Watson检验。H0: 误差项是相互独立的; H1: 误差项是相关的检验统计量：DW接近于0，表示残差中存在正自相关；如果DW接近于4，表示残差中存在负自相关；如果DW接近于2，表示残差独立性。二、R语言实现还是用lm()函数实现，不同是需要设置更复杂的formula格式：yx1+x2只考虑自变量的主效应（y=k1x1+k2x2），y.表示全部自变量的主效应；yx1+x2+x1:x2考虑主效应和交互效应（y=k1x1+k2x2+k3x1x2）；yx1*x2考虑全部主效应和交互效应的简写（效果同上）；y(x1+x2+x3)2考虑主效应以及至2阶以下的交互效应，相当于x1+x2+x3+x1:x2+x2:x3+x1:x3yx1%in%x2x1含于x2，相当于x2+x2:x1y(x1+x2)2-x1:x2表示从(x1+x2)2中去掉x1:x2yx1+I(x2+x3)2)使用I()函数，相当于用(x2+x3)2计算出新变量h，然后yx1+h function在表达式中使用数学函数，例如log(y)x1+x2三、实例 例2 现有19902009年财政收入的数据revenue.txt：各变量分别表示：y：财政收入（亿元）x1：第一产业国内生产总值（亿元）x2：第二产业国内生产总值（亿元）x3：第三产业国内生产总值（亿元）x4：人口数（万人）x5：社会消费品零售总额（亿元）x6：受灾面积（万公顷）做多元线性回归分析。setwd("E:/办公资料/R语言/R语言学习系列/codes")revenue=read.table("revenue.txt",header=TRUE) lm.reg=lm(yx1+x2+x3+x4+x5+x6,revenue)summary(lm.reg)Residuals: Min 1Q Median 3Q Max -295.71 -173.52 26.59 90.16 370.01 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.046e+04 3.211e+03 18.829 8.12e-11 *x1 -1.171e-01 8.638e-02 -1.356 0.19828 x2 3.427e-02 3.322e-02 1.032 0.32107 x3 6.182e-01 4.103e-02 15.067 1.31e-09 *x4 -5.152e-01 2.930e-02 -17.585 1.91e-10 *x5 -1.104e-01 2.878e-02 -3.837 0.00206 * x6 -1.864e-02 1.023e-02 -1.823 0.09143 . -Signif.codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Residual standard error: 234.8 on 13 degrees of freedomMultiple R-squared: 0.9999,Adjusted R-squared: 0.9999 F-statistic: 2.294e+04 on 6 and 13 DF, p-value: < 2.2e-16 说明：拟合优度R2=0.9999, 效果非常好。但是多元回归时，自变量个数越多，拟合优度必然越好，所以还要看检验效果和回归系数是否显著。结果解释、回归方程、回归预测与前文类似（略）。结合显著性代码可看出：x1和x2不显著，x6只在0.1显著水平下显著，故应考虑剔除x1和x2.R语言中提供了update()函数，用来在原模型的基础上进行修正，还可以对变量进行运算，其基本格式为：update(object, formula., ., evaluate=TRUE)其中，object为前面拟合好的原模型对象；formula指定模型的格式，原模型不变的部分用“.”表示，只写出需要修正的地方即可，例如update(lm.reg, .+x7) 表示添加一个新的变量update(lm.reg, sqrt(.).) 表示对因变量y开方，再重新拟合回归模型lm.reg2<-update(lm.reg, .-x1-x2) #剔除自变量x1,x2summary(lm.reg2)Residuals: Min 1Q Median 3Q Max -325.62 -147.54 14.07 108.28 427.42 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.339e+04 2.346e+03 27.020 3.89e-14 *x3 6.584e-01 1.548e-02 42.523 < 2e-16 *x4 -5.438e-01 1.981e-02 -27.445 3.09e-14 *x5 -1.392e-01 1.918e-02 -7.256 2.80e-06 *x6 -1.803e-02 9.788e-03 -1.842 0.0854 .-Signif.codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Residual standard error: 233.6 on 15 degrees of freedomMultiple R-squared: 0.9999,Adjusted R-squared: 0.9999 F-statistic: 3.476e+04 on 4 and 15 DF, p-value: < 2.2e-16（三）逐步回归多元线性回归模型中，并不是所有的自变量都与因变量有显著关系，有时有些自变量的作用可以忽略。这就需要考虑怎样从所有可能有关的自变量中挑选出对因变量有显著影响的部分自变量。逐步回归的基本思想是，将变量一个一个地引入或剔出，引入或剔出变量的条件是“偏相关系数”经检验是显著的，同时每引入或剔出一个变量后，对已选入模型的变量要进行逐个检验，将不显著变量剔除或将显著的变量引入，这样保证最后选入的所有自变量都是显著的。逐步回归每一步只有一个变量引入或从当前的回归模型中剔除，当没有回归因子能够引入或剔出模型时，该过程停止。R语言中，用step()函数进行逐步回归，以AIC信息准则作为选入和剔除变量的判别条件。AIC是日本统计学家赤池弘次，在熵概念的基础上建立的：AIC=2(p+1)-2ln(L)其中，p为回归模型的自变量个数，L是似然函数。注：AIC值越小越被优先选入。基本格式：step(object, direction=,steps=, k=2, .)其中，object是线性模型或广义线性模型的返回结果；direction确定逐步回归的方法，默认“both”综合向前向后法，“backward”向后法（先把全部自变量加入模型，若无统计学意义则剔出模型），“forward”向前法（先将部分自变量加入模型，再逐个添加其它自变量，若有统计学意义则选入模型）；steps表示回归的最大步数，默认1000；k默认=2, 输出为AIC值，= log(n)有时输出BIC或SBC值。 另外，有时还需要借助使用drop1(object)和add1(object)函数，其中object为逐步回归的返回结果，判断剔除或选入一个自变量，AIC值的变化情况，以筛选选入模型的自变量。lm.step<-step(lm.reg)summary(lm.step)Call:lm(formula = y x3 + x4 + x5 + x6, data = revenue)Residuals: Min 1Q Median 3Q Max -325.62 -147.54 14.07 108.28 427.42 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.339e+04 2.346e+03 27.020 3.89e-14 *x3 6.584e-01 1.548e-02 42.523 < 2e-16 *x4 -5.438e-01 1.981e-02 -27.445 3.09e-14 *x5 -1.392e-01 1.918e-02 -7.256 2.80e-06 *x6 -1.803e-02 9.788e-03 -1.842 0.0854 .-Signif.codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Residual standard error: 233.6 on 15 degrees of freedomMultiple R-squared: 0.9999,Adjusted R-squared: 0.9999 F-statistic: 3.476e+04 on 4 and 15 DF, p-value: < 2.2e-16说明：默认输出每步的结果（略），进行了3步回归，逐步剔除了自变量x1和x2，AIC值逐步减小，最终得到最优的模型。drop1(lm.step)Single term deletionsModel:y x3 + x4 + x5 + x6 Df Sum of Sq RSS AIC<none> 222.40x3 1 316.40x4 1 299.12x5 1 250.52x6 1 224.47lm.reg3<-lm(yx3+x4+x5,revenue)summary(lm.reg3)Call:lm(formula = y x3 + x4 + x5, data = revenue)Residuals: Min 1Q Median 3Q Max -336.34 -186.82 1.52 89.46 437.84 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.284e+04 2.494e+03 25.191 2.66e-14 *x3 6.614e-01 1.651e-02 40.066 < 2e-16 *x4 -5.467e-01 2.118e-02 -25.813 1.81e-14 *x5 -1.412e-01 2.053e-02 -6.877 3.72e-06 *-Signif.codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Residual standard error: 250.5 on 16 degrees of freedomMultiple R-squared: 0.9999,Adjusted R-squared: 0.9998 F-statistic: 4.032e+04 on 3 and 16 DF, p-value: < 2.2e-16 说明：使用drop1()函数考察分别剔除每个自变量，AIC值变化的情况，可以看出不剔除x6与剔除x6，AIC值只从222.40变大到224.47，相对其它自变量变化很小。所以，可以考虑剔除掉x6，重新做多元线性回归。（四）回归诊断回归分析之后，还需要从残差的随机性、强影响分析、共线性方面进行诊断。一、残差诊断(1) 残差y.res<-lm.reg3$residuals #回归模型的残差y.fit<-predict(lm.reg3) #回归模型的预测值plot(y.resy.fit,main="残差图") #绘制残差图，以预测值作为横坐标 说明：从图形看，残差分布比较均匀，大致满足随机性。shapiro.test(y.res) #残差的正态性检验Shapiro-Wilk normality testdata: y.resW = 0.94206, p-value = 0.2622说明：p值=0.2622>0.05, 接受原假设，即残差服从正态分布。(2) 标准化残差残差与数据的数量级有关，除以标准误差后得到标准化残差。理想的标准化残差服从N(0,1).rs<-rstandard(lm.reg3) #得到标准化残差plot(rsy.fit,main="标准残差图")shapiro.test(rs) #标准化残差的正态性检验Shapiro-Wilk normality testdata: rsW = 0.97766, p-value = 0.9004(3) 学生化残差为了回避标准化残差的方差齐性假设，使用学生化残差。rst<-rstudent(lm.reg3)plot(rsy.fit,main="学生化残差图")shapiro.test(rst)Shapiro-Wilk normality testdata: rstW = 0.97463, p-value = 0.848(4) 残差自相关性的Durbin-Watson检验使用car包中的函数：durbinWatsonTest(model,alternative=c("two.side", "positive", "negative")H0：序列不存在自相关性library(car)durbinWatsonTest(lm.reg3) lag Autocorrelation D-W Statistic p-value 1 -0. 2.42579 0.77 Alternative hypothesis: rho != 0二、强影响分析对参数估计或预测值有异常影响的数据，称为强影响数据。回归模型应当具有一定的稳定性，若个别一两组数据对估计有异常大的影响，剔除后将得到与原来差异很大的回归方程，从而有理由怀疑原回归方程是否真正描述了变量间的客观存在的关系。1. 反映这种强影响的统计量有4种及函数：Leveragehatvalues(model)DEFITSdffits(model)Cooks距离cooks.distance(model) COVRATIOcovratio(model)另外，influence.measures(model)函数，可以汇总上述4种统计量，判断强影响点。influence.measures(lm.reg3) Influence measures of lm(formula = y x3 + x4 + x5, data = revenue) : dfb.1_ dfb.x3 dfb.x4 dfb.x5 dffit cov.r cook.d hat inf18 0. -3.04124 -0. 2. -3.40945 0.810 2.14e+00 0.6347 *19 0. -0.09558 -0. 0. 1.61704 0.515 5.04e-01 0.3127 *20 -1. 1.56506 1. -1. -3.33452 1.426 2.25e+00 0.6996 * 说明：判断出第18, 19, 20个样本是强影响点。2. Bonferroni离群点检验使用car包中的函数outlierTest(model)library(car)outlierTest(lm.reg3)No Studentized residuals with Bonferonni p < 0.05Largest |rstudent|: rstudent unadjusted p-value Bonferonni p18 -2. 0.02064 0.4128注：去掉强影响点，重新做多元线性回归（略）。三、共线性诊断回归分析中很容易发生模型中两个或两个以上的自变量高度相关，从而引起最小二乘估计可能很不精确（称为共线性问题）。在实际中最常见的问题是一些重要的自变量很可能由于在假设检验中t值不显著而被不恰当地剔除了。共线性诊断问题就是要找出哪些变量间存在共线性关系。 1. 模型条件数检验使用函数kappa(z, exact=FALSE, )，其中，z为矩阵XTX，或lm、glm的返回对象；exact设置是否精确计算。 一般认为：当K<100时不存在多重共线性；当100K<1000时存在较强的多重共线性；当K1000时存在严重的多重共线性。x<-scale(revenue,3:8) #取出自变量数据，做标准化xx=crossprod(x) #求xx即矩阵的叉积kappa(xx)1 6132.142 2. 方差膨胀因子（VIF）检验使用car包中的函数vif(model)，该函数还能判断哪些自变量间存在共线性。一般认为：当vif<10时不存在多重共线性；当10vif<100时，存在较强的多重共线性；当vif100时存在严重的多重共线性。lm.reg<-lm(yx1+x2+x3+x4+x5+x6,revenue)vif(lm.reg) x1 x2 x3 x4 x5 x6 196. 777. 1014. 10. 342. 1. cor(revenue$x2, revenue$x3) #x2和x3的vif值最大，考察二者的相关性1 0. 可见，x2和x3存在严重的共线性，应该考虑剔除其中的一个。（五）岭回归多元线性回归分析中，我们会在众多变量中选择对因变量显著性影响大的那些自变量。但常常会遇到一个问题：在某些情况下，增加或剔除一个自变量后，回归系数变化很大甚至改变符号。主要原因就是变量之间存在多重共线性。岭回归分析是一种专用于共线性数据分析的有偏估计回归方法，实质上是一种改良的最小二乘估计法，通过放弃最小二乘法的无偏性，以损失部分信息、降低精度为代价，获得回归系数更为符合实际、更可靠的回归方法。基本原理：当自变量间存在多重共线性时，有，考虑加上一个正常数矩阵kI, (k>0), 则接近奇异的程度就会比小很多，从而消除了多重共线性。考虑到变量的量纲，应先对数据进行标准化，仍记为X，称为的岭回归估计，k称为岭参数。显然岭回归估计的值比最小二乘估计值稳定，当k=0时，岭回归就是普通的最小二乘估计。由于岭参数k不是唯一确定的，故岭回归估计实际上是回归参数的一个估计族（关于k的函数），该函数曲线称为岭迹。根据岭迹图可以选择合适的k值，称为岭迹法。其一般原则是： (1) 各回归系数的岭估计基本稳定；(2) 最小二乘估计的回归系数符合不合理时，岭估计参数的符合变的合理；(3) 回归系数每月不合乎实际意义的绝对值；(4) 残差平方和增加不太多。岭回归用MASS包中的函数lm.ridge()实现，基本格式为：lm.ridge(formula, data, subset, na.action, lambda=0,.)其中，formula为回归公式的形式，data为数据框，subset设置其子集；lambda为岭参数的标量或矢量。例3 商品销售量y，考虑4个自变量：可支配收入x1, 商品价格指数x2, 商品的社会保有量x3, 其它消费品平均价格指数x4