创业维艰，小本经营的我们也充满了对财富的渴望，绝不是奢求一夜暴富，一本万利，腰缠万贯，我们期望的经营有道，旱涝保收，恰如孟子对曰：“仰足以事父母，俯足以畜妻子，乐岁终身饱，凶年免于死亡” 。
小白狼想要的是什么？
小白狼（化名），一个城市的小个体户，他放弃了自己年薪20万的工作，每年投入40万，经营这一家餐馆 [舌尖上的中餐厅]，小白狼的期望餐馆的营业额至少要高于：60，这样自己的生活水平才不至于降低，小白狼每天焦头烂额，注意力分散在：服务人员的协调，食材的供应，薪水的发放等等，小白狼身心疲惫，觉得自己比附近的对手们做的都要努力，为什么做不好呢？当我们把注意力都打散在零散的人和事，常常会忽略自己想要的，重要的。不管你做了多少事那是都只是过程，时刻强化自己的想要的是什么。经营一家餐馆小白狼（化名）最想要的是：营业额，是不断增长的营业额。
小白狼（化名）餐厅营业额怎么样呢？【遭遇天花板】
小白狼披星戴月，大放狼性，自己辛苦一年，1到12月的的营业额（万元）：30,29,40,31,50,29,38,49,28,33,42,50；如R分析工具的计算和图表反馈，餐厅的月平均营业额为：37.42，大部分月份的营业额少于：35.50，最大值为：50.00，最小值为：28.00，我们很容易发现营业额的天花板。不管小白狼多么卖力，天花板永远是天花板。
代码语言：javascript复制> math <- c(30,29,40,31,50,29,38,49,28,33,42,50)> barplot(math)> barplot(math, col=c(1:13))
> summary(math) Min. 1st Qu. Median Mean 3rd Qu. Max. 28.00 29.75 35.50 37.42 43.75 50.00 
小白狼失落了，自己如此的努力营业额只在35.50的均值线上，上线波动，
最高的营业额才为：50，看来小白狼的餐馆经营失败了，远远低于自己预期的60万，收入缩减了一半。
找一个高的天花板
小白狼要的是60万的标准，却在50的天花板下苦命挣扎，你也许认为这是命运，但是我们为什们不能在开始之前，选择能够安放自己期望的天花板呢？我认为可以找到，也能够找到。
什么决定了天花板的高度呢？
我们尝试一个事物的变化是多种因素综合作用的结果，哪些是哪些重要的因素影响了餐厅营业额的高低呢？小白狼，彻夜长思，恍然大悟：一个餐馆的营业额（y）主要受到：周边居民人数（x1 /万人），用餐平均支出（x2 / 元/人），周边居民月平均收入（x3/ 元），周边餐馆数（x4/ 个），距市中心距离（x5/ km）等五个主要因素的影响。
怎么计算天花板的高度（R语言，多元线性回归）：
确定餐馆营业额y 与 五个x变量的线性关系，并建立线性关系模型
1,分析的数据预览：
代码语言：javascript复制 index y x1 x2 x3 x4 x5
1 1 53.2 163.0 168.6 6004 5 6.5
2 2 18.5 14.5 22.5 209 11 16.0
3 3 11.3 88.2 109.4 1919 10 18.2
4 4 84.7 151.6 277.0 7287 7 10.0
#只显示部分数据 2，假定线性模型：
2.1，绘制多个变量间的相关性图：
代码语言：javascript复制> library(corrgram)
> corrgram(example10_1[2:7], order=T, lower.panel=panel.shade,upper.panel=panel.pie,text=panel.txt)
6个变量之间的相关图（相关系数矩阵图）
解读相关系数矩阵图，请参考搜索，获取了解。蓝色表示正相关，红色表示负相关，对应颜色的饼图表示相关的度。
2.2 建立回归模型（检测报告，模型进行估计和检验用到了如下检测结果）：
检测报告：
代码语言：javascript复制#回归模型拟合
> fit1 <- lm(y~x1 x2 x3 x4 x5, data=example10_1)
> summary(fit1)
Call:
lm(formula = y ~ x1 x2 x3 x4 x5, data = example10_1)
Residuals:
Min 1Q Median 3Q Max
-16.7204 -6.0600 0.7152 3.2144 21.4805
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.2604768 10.4679833 0.407 0.68856
x1 0.1273254 0.0959790 1.327 0.20037
x2 0.1605660 0.0556834 2.884 0.00952 **
x3 0.0007636 0.0013556 0.563 0.57982
x4 -0.3331990 0.3986248 -0.836 0.41362
x5 -0.5746462 0.3087506 -1.861 0.07826 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 10.65 on 19 degrees of freedom
Multiple R-squared: 0.8518, Adjusted R-squared: 0.8128
F-statistic: 21.84 on 5 and 19 DF, p-value: 2.835e-07
#计算回归系数的置信区间
> confint(fit1, level=0.95)
2.5 % 97.5 %
(Intercept) -17.649264072 26.170217667
x1 -0.073561002 0.328211809
x2 0.044019355 0.277112598
x3 -0.002073719 0.003600932
x4 -1.167530271 0.501132297
x5 -1.220868586 0.071576251
#输出方差分析表
> anova(fit1)
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
x1 1 10508.9 10508.9 92.7389 9.625e-09 ***
x2 1 1347.1 1347.1 11.8878 0.002696 **
x3 1 85.4 85.4 0.7539 0.396074
x4 1 40.5 40.5 0.3573 0.557082
x5 1 392.5 392.5 3.4641 0.078262 .
Residuals 19 2153.0 113.3
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
从检验结果可以看出 x1，x2，x5 对 y 的总误差平方和贡献显著。2.3 回归模型方程式：
3，对模型进行估计和检验
3.1 拟合优度检验
多重决定系数是多元线性回归中回归平方和SSR 占 总平方和SST的比例，计算公式为：
它表示因变量y的总误差中被多少个自变量共同解释的比例。
为避免增加自变量而高估多重决定系数，统计学家使用样本量n和自变量的个数k 去调整 多重决定系数：
计算知：多重决定系数为：0.8518，说明日均营业额时与周边居民人数，用餐平均支出，周边居民月平均收入，周边餐馆数和距离市中心这5个自变量模型的拟合度较高。
Multiple R-squared: 0.8518, Adjusted R-squared: 0.8128
3.2估计标准误：
标准误 就是指 残差的标准差，计算公式：
计算知：估计的标准误差为：10.65，根据建立的多元线性回归方程，周边居民人数，用餐平均支出，周边居民月平均收入，周边餐馆数和距离市中心这5个自变量预测日均营业额时，平均的预测误差为10.65万元 Residual standard error: 10.65 on 19 degrees of freedom
3.3 模型的显著性检验：
包含如下检验：线性关系检验，回归系数检验。请参考 （2.2处：检测报告）诠释如下两种假设，此次省略具体说明。
3.4模型诊断：
绘制残差图诊断模型
代码语言：javascript复制> par(mfrow=c(1,2), mai=c(0.8,0.8,0.4,0.1),cex=0.8,cex.main=0.7)
> plot(fit1,which=1)
> plot(fit1,which=2)
残差图
在图中发现点2，4，16具有较大残差，残差的正态性存在问题（可以考虑重建模型，或者剔除较大残差值），如下我们剔除点 2，4 重新建立回归模型：
代码语言：javascript复制#再次绘制残差图
> par(mfrow=c(1,2), mai=c(0.8,0.8,0.4,0.1),cex=0.8,cex.main=0.7)
> plot(newfit1,which=1)
> plot(newfit1,which=2)
去掉点2，4的残差图
4，判别模型中是否存在多重共线性，处理方法
4.1 变量之间的高度相关性，造成结果混乱，影响参数正负号
4.2 识别共性和处理
代码语言：javascript复制#1，自变量之间的相关系数及其检验
> library(psych)
> corr.test(example10_1[3:7], use="complete")
Call:corr.test(x = example10_1[3:7], use = "complete")
Correlation matrix
x1 x2 x3 x4 x5
x1 1.00 0.74 0.88 -0.62 -0.28
x2 0.74 1.00 0.55 -0.54 -0.32
x3 0.88 0.55 1.00 -0.52 -0.29
x4 -0.62 -0.54 -0.52 1.00 0.10
x5 -0.28 -0.32 -0.29 0.10 1.00
Sample Size
[1] 25
Probability values (Entries above the diagonal are adjusted for multiple tests.)
x1 x2 x3 x4 x5
x1 0.00 0.00 0.00 0.01 0.47
x2 0.00 0.00 0.03 0.03 0.46
x3 0.00 0.00 0.00 0.04 0.47
x4 0.00 0.01 0.01 0.00 0.65
x5 0.18 0.12 0.16 0.65 0.00
#2，考察回归系数的显著性
#3，分析回归系数的正负号
#4，用容忍度和方差膨胀因子（VIF）,VIF 大于10 存在严重的共线性:
> library(carData)
> vif(fit1)
x1 x2 x3 x4 x5
8.233159 2.629940 5.184365 1.702361 1.174053
> 1/vif(fit1)
x1 x2 x3 x4 x5
0.1214601 0.3802368 0.1928877 0.5874195 0.8517500 如上分析数据可知，不存在严重的共线性问题
4.3 变量选择与逐步回归
建立模型之前就有选择的确定进入模型的自变量，也可以避免多重共线性问题，变量的选择方法主要有 向前选择，向后剔除，逐步回归 等。
逐步回归以 赤池信息准则AIC为选择标准，选择AIC最小的变量建立模型，计算公式为：
式中，n为样本量；p为模型中参数的个数（包括常数项）
代码语言：javascript复制#对模型fit1进行逐步回归
> fit2 <-step(fit1)
Start: AIC=123.39
y ~ x1 x2 x3 x4 x5
Df Sum of Sq RSS AIC
- x3 1 35.96 2189.0 121.81
- x4 1 79.17 2232.2 122.30
<none> 2153.0 123.39
- x1 1 199.42 2352.4 123.61
- x5 1 392.54 2545.6 125.58
- x2 1 942.22 3095.2 130.47
Step: AIC=121.81
y ~ x1 x2 x4 x5
Df Sum of Sq RSS AIC
- x4 1 78.22 2267.2 120.69
<none> 2189.0 121.81
- x5 1 445.69 2634.7 124.44
- x2 1 925.88 3114.9 128.63
- x1 1 1133.27 3322.3 130.24
Step: AIC=120.69
y ~ x1 x2 x5
Df Sum of Sq RSS AIC
<none> 2267.2 120.69
- x5 1 404.28 2671.5 122.79
- x2 1 1050.90 3318.1 128.21
- x1 1 1661.83 3929.0 132.43
>
> fit2 <- lm(y~x1 x2 x5, data=example10_1)
> summary(fit2)
Call:
lm(formula = y ~ x1 x2 x5, data = example10_1)
Residuals:
Min 1Q Median 3Q Max
-14.027 -5.361 -1.560 2.304 23.001
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.68928 6.25242 -0.270 0.78966
x1 0.19022 0.04848 3.923 0.00078 ***
x2 0.15763 0.05052 3.120 0.00518 **
x5 -0.56979 0.29445 -1.935 0.06656 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 10.39 on 21 degrees of freedom
Multiple R-squared: 0.8439, Adjusted R-squared: 0.8216
F-statistic: 37.85 on 3 and 21 DF, p-value: 1.187e-08
#逐步回归的方差分析
> anova(fit2)
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
x1 1 10508.9 10508.9 97.3392 2.452e-09 ***
x2 1 1347.1 1347.1 12.4775 0.001976 **
x5 1 404.3 404.3 3.7447 0.066558 .
Residuals 21 2267.2 108.0
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1根据R的逐步回归结果，得到最终的估计方程：
> par(mfcol=c(1,2),cex=0.7)
> plot(fit2, which=1)
> plot(fit2,which=2)
逐步回归模型的诊断图
4.4多模型的比较
完全模型，简化模型
代码语言：javascript复制#模型方差对比
> anova(fit1, fit2)
Analysis of Variance Table
Model 1: y ~ x1 x2 x3 x4 x5
Model 2: y ~ x1 x2 x5
Res.Df RSS Df Sum of Sq F Pr(>F)
1 19 2153.0
2 21 2267.2 -2 -114.17 0.5038 0.6121
#机器的交叉验证
#五折交叉验证
#AIC对比
> AIC(fit1, fit2)
df AIC
fit1 7 196.3408
fit2 5 193.6325对比数据，得知 fit2 模型的拟合的更好。
5，利用回归方程进行预测
基于多点的点估计，求出区间估计（均值的区间估计，个别值的预测区间）
代码语言：javascript复制#计算置信区间和预测区间
> x<- example10_1[,c(3,4,7)]
> pre<-predict(fit2)
> res<-residuals(fit2)
> zre<-rstandard(fit2)
> con_int<-predict(fit2,x,interval="confidence",level=0.95)
> pre_int<-predict(fit2,x,interval="prediction",level=0.95)
> mysummary<-data.frame(日均营业额=example10_1$y, 点预测值=pre,残差=res,标准化残差=zre,
置信下限=con_int[,2], 置信上限=con_int[,3], 预测下限=pre_int[,2], 预测上限=pre_int[,3])
> round(mysummary,3)
日均营业额 点预测值 残差 标准化残差 置信下限 置信上限 预测下限 预测上限
1 53.2 52.189 1.011 0.102 45.457 58.921 29.557 74.822
2 18.5 -4.501 23.001 2.359 -11.969 2.967 -27.363 18.361
3 11.3 21.963 -10.663 -1.072 15.685 28.240 -0.539 44.464
4 84.7 65.114 19.586 2.766 49.300 80.928 38.337 91.891
5 7.3 6.128 1.172 0.122 -2.283 14.540 -17.059 29.316
6 17.9 22.408 -4.508 -0.466 14.581 30.235 -0.574 45.390
7 2.5 1.278 1.222 0.125 -6.113 8.670 -21.559 24.116
8 27.3 34.719 -7.419 -0.747 28.400 41.038 12.206 57.232
9 5.9 10.628 -4.728 -0.472 4.904 16.351 -11.726 32.981
10 23.9 37.843 -13.943 -1.390 32.193 43.493 15.509 60.178
11 69.4 62.852 6.548 0.709 52.939 72.766 39.079 86.626
12 20.6 18.296 2.304 0.229 13.036 23.556 -3.943 40.535
13 1.9 -5.510 7.410 0.771 -13.694 2.674 -28.616 17.596
14 3.0 14.956 -11.956 -1.202 8.687 21.224 -7.543 37.455
15 7.3 8.860 -1.560 -0.169 -1.050 18.770 -14.913 32.632
16 46.2 29.831 16.369 1.699 21.745 37.917 6.759 52.903
17 78.8 78.016 0.784 0.112 62.105 93.927 51.182 104.850
18 11.1 13.167 -2.067 -0.209 6.420 19.915 -9.470 35.805
19 8.6 15.847 -7.247 -0.815 4.670 27.024 -8.481 40.175
20 48.9 50.727 -1.827 -0.184 44.205 57.250 28.156 73.299
21 22.1 27.461 -5.361 -0.536 21.671 33.251 5.090 49.831
22 11.1 0.233 10.867 1.354 -13.486 13.953 -25.362 25.829
23 8.6 22.627 -14.027 -1.395 17.181 28.073 0.344 44.911
24 48.9 48.961 -0.061 -0.006 42.722 55.200 26.470 71.452
25 22.1 27.005 -4.905 -0.490 21.220 32.790 4.636 49.374
#求 x1=50, x2=100, x5=10 时的日均营业额的点预测值，置信区间和预测区间（新值预测）
> x0<-data.frame(x1=50,x2=100,x5=10)
> predict(fit2, x0)
1
17.88685
> predict(fit2,x0, interval="confidence", level=0.95)
fit lwr upr
1 17.88685 10.98784 24.78585
> predict(fit2,x0, interval="prediction", level=0.95)
fit lwr upr
1 17.88685 -4.795935 40.56963
> 本文到此，小白狼 根据数据，通过R计算 找到了餐厅营业额和5个变量的线性模型了，小白狼只有5个变量数据，就可以轻易的计算出餐厅营业额的天花板了，省去了自己不断尝试的辛苦和经济损失。
金钱是宝贵的，时间是宝贵的，注意力是宝贵的，在这个不断变化的时代，我们需要用数据作为自己的眼睛，抓住关键，找到计算天花板的公式，给自己的努力放在更高的天花板之下，创造出真正自己期望的价值。
