因变量是分类的，则为分类分析，因变量是连续的，称为回归分析。回归分析通过建立函数表达式，用一个或者多个自变量的变化解释或预测因变量，常用于描述、探索和检验自变量和因变量之间因果关系，根据自变量的变化预测因变量的取值。通常按照自变量的个数划分为一元回归和多元回归。按照函数表达式的形式，分为线性回归和非线性回归。

一元线性回归

假设有两个变量X和Y，X为自变量，Y为因变量。则一元线性回归模型的基本结构形式为\[Y=\beta _{0}+\beta _{1}X+\varepsilon \]。\(\varepsilon\)是误差项，表示未知或不易测量的随机因素对因变量影响的总和。\(\beta _{0}\)为回归的常数，\(\beta _{1}\)为回归系数，它表示增加和减少一个单位时，Y的平均变化量。线性回归就是根据已经观测到的样本数据，应用最小二乘法获得对\(\beta _{0}\)和\(\beta _{1}\)的估计，进而得到回归方程。由于参数估计时并不知道是否存在线性关系，回归方程在应用前需要完成对回归方程的检验，即对回归模型的系数是否为零进行检验。常用t检验、F检验和相关系数检验评价回归方程的回归系数是否为0。 线性回归应用有四个前提条件：线性、独立、正态和方差齐性。线性指自变量和因变量在散点图大致呈直线趋势。独立性值观察值之间应相互独立。正态性指残差应符合正态分布。方差齐性指在自变量取值范围内，对于自变量的取值，因变量都有相同的方差。

例1 某医生分别采用盐析法和结合法测定正常皮肤中胶原蛋白的含量。盐析法只能部分提纯，结合法较为复杂单精确。该医生欲寻求找盐析法和结合法之间的关系，以便通过盐析法预测结合法的测定值。

编号 盐析法 结合法
1 6.8 546
2 7.8 553
3 8.7 562
4 8.7 563
5 8.9 570
6 9.5 575
7 10.1 581
8 10.2 605
9 10.3 607
10 10.4 621
11 11.1 624
12 12.4 626
13 13.3 632
14 13.1 640
15 13.2 656

解 做散点图，观察是否存在线性关系。由于线性回归的条件（线性、正态性、方差齐性和独立性）是通过残差来完成的，可先建立回归方程，然后通过回归诊断来完成线性回归条件的检验。 线性条件可通过散点图直接观察。

y <- c(546,553,562,563,570,575,581,605,607,621,624,626,632,640,656)
x <- c(6.8,7.8,8.7,8.7,8.9,9.5,10.1,10.2,10.3,10.4,11.1,12.4,13.3,13.1,13.2)
df <- as.data.frame(cbind(x,y))

#ggplot(df,aes(x,y))+geom_point()+stat_smooth(method = "lm")
plot(x,y)

通过散点图，可以发现二者近似直线关系，符合线性条件。

对数据进行线性回归分析。

model <- lm(y~x) #~左边为响应变量，右边为各个预测变量，预测变量之间用+符号分隔
summary(model)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.140  -8.372  -3.824   9.429  21.942 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  426.625     15.915   26.81 9.18e-13 ***
## x             16.580      1.518   10.92 6.43e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.42 on 13 degrees of freedom
## Multiple R-squared:  0.9017, Adjusted R-squared:  0.8941 
## F-statistic: 119.2 on 1 and 13 DF,  p-value: 6.426e-08

结果中Coefficients是参数估计的结果。结果显示，截距项和回归系统均有统计学意义。模型简单评价\(R^{2}\)为0.9017，校正决定系数\(R_{adj}^{2}\)为0.9841。决定系数越大表明自变量对因变量的解释程度越高。F检验检验所有的预测变量预测响应变量是否都在某个几率水平之上，结果表明(F=119.2,P=6.426e-08),方程总体有统计学意义。残差的标准误则可认为模型用自变量预测因变量的平均误差。 所建立的方程为 Y=426.625+16.580×X 对估计值做出区间估计

confint(model)
##                 2.5 %    97.5 %
## (Intercept) 392.24255 461.00814
## x            13.29973  19.86039
plot(df$x,df$y)
abline(model)

###模型评价

回归诊断

主要包括三方面：（1）误差项是否满足独立性、等方差性和正态性，选择模型是否合适。（2）是否存在异常样本，回归分析的结果是否对某些样本依赖过重，即回归模型是否具备稳定性。（3）自变量之间是否存在高度相关，即是否有多重共线性问题存在。

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

par(mfrow=c(1,1))

标准方法 正态性 当预测变量值固定时，因变量成正态颁，则残差图也应是一个均值为0的正态颁。正态Q-Q图是在正态颁对应的值上，标准化残差的概率图，若满足正态假设，则图上的点应该落在45度角的直线上，若不是，则违反了正态性假设。第二幅Normal QQ-plot图中数据点分布趋于一条直线, 说明残差是服从正态分布的。

独立性 无法从图中分辨因变量值是否相互独立，只能从收集的数据中验证。本例中适用结合法进行测量时，无理由相信一个测量结果会影响另外一次的测量。

线性 若因变量与自变量线性相关，则残差值与预测（拟合）值就没有系统关联，若存在关系，则说明可能城要对回归模型进行调整。第一幅图Residual vs fitted为拟合值y对残差的图形, 可以看出,数据点都基本均匀地分布在直线y=0的两侧, 无明显趋势，满足线性假设。

方差齐性 若满足不变方差假设，则在第三幅图位置尺度图（Scale-Location Graph）中，水平线周围的点应随机分布，Scale-Location 图显示了标准化残差(standardized residuals)的平方根的分布情况，最高点为残差最大值点。第三副图显示基本符合方差齐性的要求。

第四幅图（Residuals vs Leverage）提供了单个观测点的信息，从图中可以鉴别离群点、高高杆值点和强影响点。

改进方法

正态性 通过对残差正态性检验予以证实。

model <- lm(y~x)
shapiro.test(residuals(model))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.92987, p-value = 0.2716

正态性检验结果表明W值为0.9298668，P值为0.2716282,残差符合正态分布。

独立性 判断因变量（或残差）最好的方法时依据收集数据的方式的先验知识。lmtest包提供了dwtest检验函数，car包提供了Durbin-Watson检验的函数，都能够检验误差序列的相关性。

dwtest(model)
## 
##  Durbin-Watson test
## 
## data:  model
## DW = 1.0041, p-value = 0.006564
## alternative hypothesis: true autocorrelation is greater than 0
#durbinWatsonTest(model)

本例结果P值比较显著，但根据先验知识，并不能否定因变量的独立性。

线性 可通过成分残差图(component plus residual plot)即偏残差图(partial residual plot)，判断因变量与自变量之间是否呈非线性关系，也可以看是否不同于已设定线性模型的系统偏差，图形可用car包中crPlots()函数绘制。图形存在非线性，则说明可能对预测变量的函数形式建模不够充分.

crPlots(model)

图形呈现线性，建模比较充分。car包中提供了一个linearHypothesis()函数可以自动的进行线性假设检验,比图形更为精准。根据对模型的设定，这个函数既可以用一般的方法或调整后的协方差矩阵进行F或Wald检验。

linearHypothesis(model, "x=0") #x的系数是否为0
## Linear hypothesis test
## 
## Hypothesis:
## x = 0
## 
## Model 1: restricted model
## Model 2: y ~ x
## 
##   Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
## 1     14 17249.6                                  
## 2     13  1695.8  1     15554 119.23 6.426e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# tests Beta1 = Beta2
#linear.hypothesis(fit,"x1 = x2") 
# Tests  Beta0 = Beta1 = Beta2= 1
#linear.hypothesis(fit,c("(Intercept)", "x1","x2"),rep(1,3)) 
# Tests  Beta0 = Beta1 = Beta2 = 0
#linear.hypothesis(fit,c("(Intercept)", "x1","x2"),rep(0,3)) 
# Tests Beta1 = Beta2 = 0
#linear.hypothesis(fit,c("x1","x2"),rep(0,2))

P值小于0.05，可以认为x的系数不为0。

方差齐性 通过以因变量为x轴，学生化残差为y轴做残差图，进行判断。

plot(predict(model),rstudent(model))

所有的学生化残差均在\(\pm 2\)在范围内波动，没有明显的上升或下降趋势，可以认为符合方差齐性。还通过自变量与残差绝对值的等级相关检验来判断。

spearman.test(x,abs(residuals(model))) #R中pspearman包中的spearman.
## Warning in spearman.test(x, abs(residuals(model))): Cannot compute exact p-
## values with ties
## 
##  Spearman's rank correlation rho
## 
## data:  x and abs(residuals(model))
## S = 294.76, p-value = 0.0768
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4736374
#test函数可以完成斯皮尔曼等级相关检验

#cor.test(x,abs(residuals(model)),method="spearman") #或者用cor.test()

自变量与残差绝对值的等级相关系数为0.0745175,P值大于0.05，无统计学意义，可以认为残差方差齐性。

car包提供了两个有用的函数，可判断误差方差是否恒定，ncvTest()函数生成一个计分检验，零假设为误差方差不变,备择假设为误差方差随着拟合值水平的变化而变化。若检验显著，择说明存在异方差性。spreadLevelPlot()函数创建一个添加了最佳拟合曲线的散点图，展示标准化残差绝对值与拟合值的关系。

ncvTest(model)  
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 0.4273753    Df = 1     p = 0.5132797
spreadLevelPlot(model) 

## 
## Suggested power transformation:  -2.523794

计分检验结果不显著，说明满足方差齐性的假设。通过水平分布图，可以看到其中的点在水平的最佳拟合曲线周围呈水平随机分布。若违反了该假设，将会呈现一个非水平的曲线。

线性模型假设的综合验证

gvlma(model)
## 
## Call:
## lm(formula = y ~ x)
## 
## Coefficients:
## (Intercept)            x  
##      426.63        16.58  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = model) 
## 
##                     Value p-value                Decision
## Global Stat        2.1559  0.7071 Assumptions acceptable.
## Skewness           0.4939  0.4822 Assumptions acceptable.
## Kurtosis           0.5558  0.4560 Assumptions acceptable.
## Link Function      0.3285  0.5665 Assumptions acceptable.
## Heteroscedasticity 0.7777  0.3778 Assumptions acceptable.

Global Stat给模型假设提供了一个单独的综合检验（通过/不通过），本例中，可以看到数据满足线性回归模型所有的统计假设（P=0.7071）。同时还对偏度(Skewness)、峰度(Kurtosis)、连接函数(Link function)和异方差性(Heteroscedasticity)进行了评价。

影响分析

是探查对估计有异常影响的数据，如果一个样本不服从某个模型，其余数据服从这个模型，则称该样本点为强影响点。影响分析就是区分这样的样本数据。

离群点

指那些模型预测效果不佳的观测点，通常有很大的、或正或负的残差，正残差说明模型低估了响应值，负残差说明高佑了响应值。

outlierTest(model) #Bonferroni离群点检验
## 
## No Studentized residuals with Bonferonni p < 0.05
## Largest |rstudent|:
##    rstudent unadjusted p-value Bonferonni p
## 10 2.290689           0.040881      0.61321
qqPlot(model,labels=row.names(df),id.method = "identify",simulate=T,main="QQPlot") #car包

outlierTest（）函数是根据单个最大（或正或负）残差值的显著性来判断是否有离群点，若不显著，则说明数据集中没有离群点，若显著，则必须删除该离群点，然后再检验是否还有其他离群点存在。qqPlot图中落在置信区间带外的点可被认为时离群点。本例中未发现有离群点。

高杠杆值点

是与其他预测变量有关的离群点，即它们是由许多异常的预测变量组合起来的，与响应变量值没有关系。 高杠杆值的观测点可通过帽子矩阵的值（hat statistic）判断。对于一个给定的数据集，帽子均值为p/n，其中p是模型估计的参数数目（包含截距项），n是样本量。一般来说，若观测点的帽子值大于帽子均值的2或3倍，则可认定为高杠杆值点。

 hat.plot<-function(fit){  
  p<-length(coefficients(fit))  
  n<-length(fitted(fit))  
  plot(hatvalues(fit),main="Index Plot of Hat Values")  
  abline(h=c(2,3)*p/n,col="red",lty=2)  
  identify(1:n,hatvalues(fit),names(hatvalues(fit)))  
}  
hat.plot(model) 

## integer(0)

此图中可以看到1号点是高杠杆值点。

强影响点

强影响点，即对模型参数估计值影响有些比例失衡的点。例如，当移除 模型的一个观测点时模型会发生巨大的改变，那么需要检测一下数据中是否存在强影响点。Cook距离，或称为D统计量。Cook’s D值大于4/(n-k-1)，则表明它是强影响点，其中n为样本量大小，k是预测变量数目（有助于鉴别强影响点，但并不提供关于这些点如何影响模型的信息）。

plot(model,which=4)

Cook距离(Cook’s distance)图显示了对回归的影响点。根据Cook距离，13号点可能是个强影响点。

帽子统计量、DFFITS准测、Cook统计量和COVRATIO准则在R软件可分别通过hatvalues(),dffits(),cooks.distance()和covration()函数计算。influence.measures()可对一次获得这四个统计量的结果。 影响分析综合分析

influencePlot(model) 

##       StudRes        Hat     CookD
## 1   0.6709919 0.28317427 0.3047231
## 10  2.2906887 0.06684341 0.3763698
## 13 -1.5931097 0.22573347 0.5751787
#car包中的influencePlot（）函数，可将离群点、
#杠杆点和强影响点的信息整合到一幅图形中
influence.measures(model)
## Influence measures of
##   lm(formula = y ~ x) :
## 
##     dfb.1_     dfb.x  dffit cov.r  cook.d    hat inf
## 1   0.4003 -3.69e-01  0.422 1.521 0.09286 0.2832   *
## 2  -0.1133  1.01e-01 -0.127 1.409 0.00872 0.1771    
## 3  -0.2217  1.84e-01 -0.289 1.187 0.04281 0.1119    
## 4  -0.1956  1.62e-01 -0.255 1.215 0.03370 0.1119    
## 5  -0.0910  7.34e-02 -0.125 1.276 0.00843 0.1013    
## 6  -0.1305  9.11e-02 -0.239 1.141 0.02934 0.0780    
## 7  -0.0924  3.32e-02 -0.324 1.001 0.05083 0.0674    
## 8   0.0523 -1.14e-02  0.222 1.125 0.02522 0.0668    
## 9   0.0427  1.44e-16  0.230 1.115 0.02703 0.0667    
## 10  0.0825  3.15e-02  0.613 0.609 0.14165 0.0668    
## 11 -0.0732  1.37e-01  0.361 1.000 0.06253 0.0780    
## 12  0.1404 -1.73e-01 -0.236 1.300 0.02929 0.1446    
## 13  0.6230 -7.22e-01 -0.860 1.033 0.33083 0.2257    
## 14  0.1294 -1.51e-01 -0.184 1.445 0.01821 0.2052    
## 15 -0.3898  4.54e-01  0.546 1.257 0.14826 0.2153

纵坐标超过2或小于-2的州可被认为是离群点，水平轴超过0.2或0.3的州有高杠杆值（通常为预测值的组合）。圆圈大小与影响成比例，圆圈很大的点可能是对模型估计造成的不成比例影响的强影响点。influence.measures()的inf用×标注异常值。

共线性，条件数

本例只有一个自变量，不涉及。

预测新值及其置信区间

把预测变量数据保存为一个数据框，调用predict函数，将数据框做为参数。

preds <- data.frame(x=14)
#默认0.95的置信水平，可通过level改变
predict(model,newdata = preds,interval = "prediction") 
##        fit      lwr      upr
## 1 658.7462 630.5198 686.9727

改进措施

model2 <- lm(y~x,subset=-1)
gvlma(model2)
## 
## Call:
## lm(formula = y ~ x, subset = -1)
## 
## Coefficients:
## (Intercept)            x  
##      420.12        17.15  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = model2) 
## 
##                     Value p-value                Decision
## Global Stat        4.4852  0.3443 Assumptions acceptable.
## Skewness           0.6001  0.4385 Assumptions acceptable.
## Kurtosis           0.3603  0.5484 Assumptions acceptable.
## Link Function      2.6142  0.1059 Assumptions acceptable.
## Heteroscedasticity 0.9106  0.3399 Assumptions acceptable.
influence.measures(model2)
## Influence measures of
##   lm(formula = y ~ x, subset = -1) :
## 
##     dfb.1_   dfb.x   dffit cov.r  cook.d    hat inf
## 2  -0.0447  0.0410 -0.0486 1.575 0.00129 0.2455   *
## 3  -0.2319  0.2028 -0.2801 1.294 0.04113 0.1502    
## 4  -0.1993  0.1743 -0.2407 1.321 0.03068 0.1502    
## 5  -0.0766  0.0659 -0.0964 1.359 0.00504 0.1341    
## 6  -0.1506  0.1193 -0.2330 1.204 0.02831 0.0968    
## 7  -0.1288  0.0790 -0.3191 1.041 0.04993 0.0761    
## 8   0.0877 -0.0483  0.2480 1.123 0.03137 0.0742    
## 9   0.0763 -0.0355  0.2529 1.113 0.03248 0.0729    
## 10  0.1620 -0.0554  0.6536 0.569 0.15525 0.0719    
## 11 -0.0490  0.1069  0.3587 0.998 0.06170 0.0784    
## 12  0.1555 -0.1875 -0.2589 1.309 0.03533 0.1502    
## 13  0.7317 -0.8307 -0.9866 0.975 0.41740 0.2455    
## 14  0.1725 -0.1977 -0.2404 1.473 0.03094 0.2211    
## 15 -0.3708  0.4228  0.5078 1.337 0.13056 0.2330

提示2号点也是是个异常值。

model3 <- lm(y~x,subset=c(-1,-2))
gvlma(model3)
## 
## Call:
## lm(formula = y ~ x, subset = c(-1, -2))
## 
## Coefficients:
## (Intercept)            x  
##      421.00        17.08  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = model3) 
## 
##                     Value p-value                   Decision
## Global Stat        5.4865 0.24092    Assumptions acceptable.
## Skewness           0.4898 0.48401    Assumptions acceptable.
## Kurtosis           0.5080 0.47599    Assumptions acceptable.
## Link Function      4.0988 0.04291 Assumptions NOT satisfied!
## Heteroscedasticity 0.3899 0.53236    Assumptions acceptable.
influence.measures(model3)
## Influence measures of
##   lm(formula = y ~ x, subset = c(-1, -2)) :
## 
##     dfb.1_   dfb.x  dffit cov.r cook.d    hat inf
## 3  -0.2902  0.2616 -0.335 1.379 0.0590 0.1973    
## 4  -0.2505  0.2258 -0.289 1.410 0.0444 0.1973    
## 5  -0.0997  0.0889 -0.119 1.447 0.0077 0.1751    
## 6  -0.1891  0.1594 -0.262 1.251 0.0361 0.1220    
## 7  -0.1722  0.1261 -0.338 1.065 0.0564 0.0893    
## 8   0.1177 -0.0821  0.255 1.158 0.0334 0.0859    
## 9   0.1062 -0.0694  0.257 1.145 0.0340 0.0830    
## 10  0.2387 -0.1422  0.664 0.583 0.1612 0.0806    
## 11 -0.0175  0.0695  0.345 1.018 0.0577 0.0802    
## 12  0.1462 -0.1748 -0.248 1.336 0.0327 0.1529    
## 13  0.7307 -0.8209 -0.979 1.007 0.4137 0.2594    
## 14  0.1677 -0.1901 -0.233 1.520 0.0292 0.2318    
## 15 -0.3759  0.4241  0.512 1.373 0.1334 0.2453

删除1、2观测值后，模型的影响分析的结果变得更好。但应该对删除观测点的方法谨慎，因为收集数据的异常点可能是最有意义东西，除非确定数据点时记录错误或者没有相关遵守规程。

多元线性回归

例２　某项“冠状动脉缓慢血流现象”的影响因素的研究，以前降支、回旋支、右冠状动脉三支血管的平均TIMI帧基数(MTFC)表示，调查的影响因素有年龄(AGE,岁)、收缩压(SBP,mmHg)、舒张压(DBP,mmHg)、白细胞(WBC,\(10^{2}\)/L),寻找影响MTFC变化的因素。

age sbp dbp wbc mtfc
43 110 50 6.19 33.67
63 105 60 6.03 26.67
59 100 60 5.28 23
78 100 60 6.52 26
67 100 60 7.31 28
65 119 61 5.67 30.33
66 120 64 5.11 27
73 130 88 6.40 47
53 113 68 4.41 27.67
76 120 70 4.20 37.33
76 136 70 5.38 35.67
76 130 70 4.94 31.33
68 126 70 4.56 32.33
61 136 70 5.42 30.67
78 124 70 5.75 37.67
80 110 70 4.68 36
74 140 70 8.67 41
75 130 70 6.62 41.67
66 130 70 6.86 22
55 114 70 7.52 23.33
71 120 70 4.94 25.67
62 130 70 4.59 25
69 130 70 4.26 27
45 110 70 10.21 29
79 120 70 6.46 30.33
58 110 70 4.70 27
65 100 70 6.06 28
44 119 70 5.55 22.33
53 110 70 14.0 29.33
62 130 72 7.29 43
62 118 72 3.97 27.33
53 122 74 3.97 18.33
71 130 75 3.78 31
54 116 75 4.35 22.33
64 120 76 6.59 30
71 140 78 5.70 35.67
50 121 78 5.27 40.33
51 138 80 5.65 34.67
73 130 80 7.45 35.33
64 138 80 6.58 33.67
40 130 80 7.51 35.33
72 120 80 4.42 34
51 100 80 7.85 21
49 120 89 5.14 20.67
63 150 90 8.18 42.67
56 130 90 5.23 30.67
69 160 90 7.10 39
78 130 90 6.03 29
78 120 90 4.52 30.67
61 150 92 7.52 40
76 142 92 4.66 38
51 140 100 5.70 28.33
51 140 100 6.71 42.67
57 160 100 6.14 41
63 190 100 5.25 46
69 150 80 6.33 22.67 
records <- read.table("example1")
ex <- rename(records, c("V1"="age","V2"="sbp","V3"="dbp","V4"="wbc","V5"="mtfc"))
attach(ex)
#scatterplotMatrix()函数默认在非对角线区域绘制变量间的散点图，并添加平滑（loess）和线性拟合区间
scatterplotMatrix(ex,spread=F,lty.smooth=2,main="scatter plot matrix")

绘制因变量与自变量的散点图矩阵显示mtfc和wbc线性关系不是很好。首先做单因素的线性回归，尽管有时候单因素的分析不是必须的，其结果也不一定可靠，但有助于初步探索自变量和因变量之间的关系。

summary(lm(mtfc~age))
## 
## Call:
## lm(formula = mtfc ~ age)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.666  -5.065  -1.004   4.454  14.382 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 21.39688    5.62809   3.802 0.000367 ***
## age          0.16225    0.08742   1.856 0.068930 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.86 on 54 degrees of freedom
## Multiple R-squared:  0.05996,    Adjusted R-squared:  0.04255 
## F-statistic: 3.444 on 1 and 54 DF,  p-value: 0.06893
summary(lm(mtfc~sbp))
## 
## Call:
## lm(formula = mtfc ~ sbp)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.6165  -3.2236   0.1771   2.7278  14.4407 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.83276    5.87395   0.312    0.756    
## sbp          0.23636    0.04607   5.130 4.04e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.802 on 54 degrees of freedom
## Multiple R-squared:  0.3277, Adjusted R-squared:  0.3152 
## F-statistic: 26.32 on 1 and 54 DF,  p-value: 4.038e-06
summary(lm(mtfc~dbp))
## 
## Call:
## lm(formula = mtfc ~ dbp)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.2937  -4.5413   0.2055   3.9111  12.2893 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 12.69827    5.88525   2.158  0.03542 * 
## dbp          0.25017    0.07663   3.265  0.00191 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.466 on 54 degrees of freedom
## Multiple R-squared:  0.1648, Adjusted R-squared:  0.1494 
## F-statistic: 10.66 on 1 and 54 DF,  p-value: 0.001906
summary(lm(mtfc~wbc))
## 
## Call:
## lm(formula = mtfc ~ wbc)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.0883  -5.1099  -0.5967   5.3872  15.0602 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  27.9326     3.4599   8.073 7.51e-11 ***
## wbc           0.6261     0.5533   1.132    0.263    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.993 on 54 degrees of freedom
## Multiple R-squared:  0.02317,    Adjusted R-squared:  0.005077 
## F-statistic: 1.281 on 1 and 54 DF,  p-value: 0.2628
summary(lm(mtfc~age+sbp+dbp+wbc))
## 
## Call:
## lm(formula = mtfc ~ age + sbp + dbp + wbc)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.3010  -3.0202  -0.0171   2.3626  12.5424 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -13.77139    8.23278  -1.673  0.10050   
## age           0.16272    0.07463   2.180  0.03388 * 
## sbp           0.20813    0.06273   3.318  0.00168 **
## dbp           0.04392    0.09327   0.471  0.63971   
## wbc           0.91334    0.45334   2.015  0.04922 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.586 on 51 degrees of freedom
## Multiple R-squared:  0.4114, Adjusted R-squared:  0.3652 
## F-statistic: 8.911 on 4 and 51 DF,  p-value: 1.552e-05

分析结果表明，单因素分析中dbp有统计学，而多因素分析中却没有统计学意义。分析自变量的相关系数

cor(ex[1:4])
##            age         sbp          dbp          wbc
## age  1.0000000  0.10515656 -0.059799697 -0.222746594
## sbp  0.1051566  1.00000000  0.691707562 -0.030466052
## dbp -0.0597997  0.69170756  1.000000000  0.003302349
## wbc -0.2227466 -0.03046605  0.003302349  1.000000000
cor.test(dbp,sbp)
## 
##  Pearson's product-moment correlation
## 
## data:  dbp and sbp
## t = 7.0384, df = 54, p-value = 3.569e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5241180 0.8077234
## sample estimates:
##       cor 
## 0.6917076

相关分析结果表明sbp和dbp有明显的正相关作用，说明单因素分析中dbp对因变量的作用同时包含了部分sbp的正向作用。在删除dbp变量，继续建模。

summary(lm(mtfc~age+sbp+wbc))
## 
## Call:
## lm(formula = mtfc ~ age + sbp + wbc)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.5790  -3.0970  -0.0079   2.2345  12.6560 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -12.64491    7.81842  -1.617   0.1119    
## age           0.15634    0.07284   2.146   0.0365 *  
## sbp           0.22890    0.04427   5.170 3.79e-06 ***
## wbc           0.91178    0.44992   2.027   0.0479 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.544 on 52 degrees of freedom
## Multiple R-squared:  0.4088, Adjusted R-squared:  0.3747 
## F-statistic: 11.99 on 3 and 52 DF,  p-value: 4.441e-06

结果显示三个因素均有统计学意义，F检验也通过了，但决定系数较低，说明自变量对因变量的解释程度较低。查看wbc变量估计结果，其标准误为0.44远高于age和sbp变量，计算这三个变量的变异系数。

cv <- function(x){
  return(100*sd(x)/mean(x))
}

cv(age)
## [1] 16.65874
cv(sbp)
## [1] 13.43599
cv(wbc)
## [1] 28.30682

wbc变量的变异系数远高于其他连个变量，为减少wbc变量的变异，对其进行对数变换后重新建模

fit <- lm(mtfc~age+sbp+log10(wbc))
summary(fit)
## 
## Call:
## lm(formula = mtfc ~ age + sbp + log10(wbc))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.7634  -3.0155   0.0769   2.4431  12.6482 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -18.67195    9.17172  -2.036   0.0469 *  
## age           0.15739    0.07214   2.182   0.0337 *  
## sbp           0.22466    0.04398   5.108 4.71e-06 ***
## log10(wbc)   15.65527    7.04993   2.221   0.0308 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.504 on 52 degrees of freedom
## Multiple R-squared:  0.4174, Adjusted R-squared:  0.3838 
## F-statistic: 12.42 on 3 and 52 DF,  p-value: 3.068e-06
gvlma(fit)
## 
## Call:
## lm(formula = mtfc ~ age + sbp + log10(wbc))
## 
## Coefficients:
## (Intercept)          age          sbp   log10(wbc)  
##    -18.6720       0.1574       0.2247      15.6553  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = fit) 
## 
##                      Value p-value                Decision
## Global Stat        2.35362  0.6710 Assumptions acceptable.
## Skewness           0.01377  0.9066 Assumptions acceptable.
## Kurtosis           1.70854  0.1912 Assumptions acceptable.
## Link Function      0.05547  0.8138 Assumptions acceptable.
## Heteroscedasticity 0.57583  0.4479 Assumptions acceptable.
influence.measures(fit)
## Influence measures of
##   lm(formula = mtfc ~ age + sbp + log10(wbc)) :
## 
##       dfb.1_   dfb.age   dfb.sbp  dfb.l10.    dffit cov.r   cook.d    hat
## 1   3.59e-01 -0.428546 -0.177294 -0.031666  0.53990 0.973 7.06e-02 0.0972
## 2  -6.61e-03 -0.001149  0.012231 -0.001698 -0.01559 1.134 6.19e-05 0.0474
## 3  -5.16e-02  0.012767  0.054344  0.016056 -0.06935 1.152 1.22e-03 0.0666
## 4   1.45e-02 -0.130260  0.134767 -0.063314 -0.19695 1.196 9.83e-03 0.1170
## 5   3.30e-04  0.002551 -0.005895  0.003831  0.00787 1.183 1.58e-05 0.0867
## 6   1.82e-03  0.000999 -0.002664 -0.000403  0.00653 1.105 1.09e-05 0.0220
## 7  -3.52e-02 -0.011663  0.026555  0.033058 -0.08254 1.088 1.73e-03 0.0260
## 8  -2.65e-01  0.321937  0.032506  0.187859  0.47948 0.735 5.27e-02 0.0387
## 9   1.22e-01 -0.078254 -0.042390 -0.088898  0.13676 1.150 4.75e-03 0.0769
## 10  9.96e-02  0.188811 -0.089872 -0.210711  0.37910 1.000 3.53e-02 0.0691
## 11 -6.63e-03  0.010612  0.004404 -0.000950  0.01578 1.134 6.34e-05 0.0469
## 12  1.30e-02 -0.054026 -0.005469  0.022504 -0.08282 1.121 1.74e-03 0.0468
## 13  2.23e-02  0.009516 -0.001399 -0.038892  0.05995 1.114 9.14e-04 0.0366
## 14 -6.64e-03  0.021360 -0.035047  0.021462 -0.07129 1.096 1.29e-03 0.0279
## 15 -6.59e-02  0.157256 -0.032776  0.027252  0.19287 1.087 9.36e-03 0.0545
## 16  5.78e-02  0.283320 -0.207803 -0.099519  0.41625 1.039 4.27e-02 0.0905
## 17 -1.03e-01  0.065269  0.032520  0.093370  0.12315 1.195 3.86e-03 0.1050
## 18 -1.69e-01  0.201864  0.012127  0.126143  0.27755 1.015 1.91e-02 0.0497
## 19  1.95e-01 -0.113143 -0.051565 -0.225404 -0.39696 0.749 3.64e-02 0.0288
## 20 -3.45e-02  0.081568  0.105501 -0.140112 -0.26248 1.035 1.71e-02 0.0521
## 21 -3.93e-02 -0.073744  0.051051  0.060842 -0.16475 1.058 6.82e-03 0.0355
## 22 -1.02e-01  0.054556 -0.038540  0.147740 -0.20919 1.032 1.09e-02 0.0383
## 23 -5.86e-02 -0.025869 -0.022187  0.129706 -0.17886 1.083 8.06e-03 0.0489
## 24  9.71e-05 -0.002775 -0.001969  0.004529  0.00672 1.279 1.15e-05 0.1553
## 25  6.80e-02 -0.132009  0.045162 -0.060857 -0.16090 1.135 6.56e-03 0.0723
## 26  5.00e-02 -0.020702 -0.029534 -0.032452  0.05844 1.138 8.70e-04 0.0550
## 27  3.22e-02  0.015698 -0.071003  0.011109  0.08441 1.146 1.81e-03 0.0645
## 28 -1.73e-01  0.218065  0.025369  0.067449 -0.25220 1.124 1.60e-02 0.0870
## 29  1.25e-01  0.015424  0.098180 -0.340615 -0.37780 1.424 3.61e-02 0.2640
## 30 -1.36e-01  0.005373  0.046873  0.210047  0.32210 0.889 2.50e-02 0.0334
## 31  1.39e-02 -0.003987 -0.003950 -0.014962  0.01822 1.160 8.46e-05 0.0688
## 32 -4.12e-01  0.295686  0.019636  0.398931 -0.50812 0.988 6.28e-02 0.0945
## 33  5.00e-03  0.002158  0.001317 -0.010801  0.01348 1.170 4.63e-05 0.0763
## 34 -1.61e-01  0.103681  0.043700  0.128725 -0.18579 1.125 8.72e-03 0.0721
## 35  2.54e-03 -0.005733  0.011809 -0.015582 -0.03493 1.105 3.11e-04 0.0257
## 36  2.07e-03 -0.001852 -0.002187 -0.000092 -0.00425 1.122 4.61e-06 0.0367
## 37  4.33e-01 -0.470455 -0.054442 -0.234302  0.60734 0.726 8.39e-02 0.0566
## 38  3.13e-02 -0.086518  0.054277 -0.026859  0.11746 1.128 3.50e-03 0.0584
## 39  1.10e-02 -0.009948 -0.000679 -0.010731 -0.01611 1.148 6.61e-05 0.0588
## 40  2.50e-02 -0.002932 -0.025733 -0.018735 -0.05022 1.108 6.42e-04 0.0307
## 41  7.12e-02 -0.277512  0.056951  0.073964  0.33680 1.145 2.84e-02 0.1170
## 42  6.45e-02  0.068825 -0.046912 -0.104893  0.18221 1.083 8.36e-03 0.0500
## 43 -9.97e-02  0.103237  0.193550 -0.135239 -0.30819 1.124 2.38e-02 0.1001
## 44 -2.40e-01  0.251372  0.034486  0.136304 -0.32085 1.030 2.55e-02 0.0648
## 45 -1.44e-01  0.008134  0.124125  0.122439  0.19957 1.144 1.01e-02 0.0862
## 46  1.25e-03 -0.001574  0.000579 -0.001129  0.00258 1.120 1.70e-06 0.0353
## 47  1.21e-01 -0.032307 -0.127150 -0.056930 -0.15983 1.186 6.48e-03 0.1041
## 48  1.43e-01 -0.221494 -0.008492 -0.067026 -0.27338 1.036 1.86e-02 0.0553
## 49 -7.84e-04 -0.004803  0.001934  0.002891 -0.00749 1.158 1.43e-05 0.0668
## 50 -5.33e-02 -0.007733  0.060126  0.041048  0.08624 1.156 1.89e-03 0.0721
## 51 -2.83e-02  0.055821  0.050276 -0.042772  0.11402 1.137 3.30e-03 0.0638
## 52 -4.39e-02  0.148687 -0.106041  0.042380 -0.20664 1.095 1.07e-02 0.0617
## 53 -1.44e-02 -0.280438  0.218036  0.069937  0.43272 0.922 4.52e-02 0.0607
## 54 -5.01e-02 -0.052872  0.134352  0.000523  0.15460 1.183 6.07e-03 0.1016
## 55 -4.31e-02 -0.015028  0.102207 -0.015656  0.10632 1.504 2.88e-03 0.2829
## 56  5.28e-01 -0.200230 -0.591615 -0.178380 -0.79877 0.551 1.35e-01 0.0585
##    inf
## 1     
## 2     
## 3     
## 4     
## 5     
## 6     
## 7     
## 8    *
## 9     
## 10    
## 11    
## 12    
## 13    
## 14    
## 15    
## 16    
## 17    
## 18    
## 19   *
## 20    
## 21    
## 22    
## 23    
## 24   *
## 25    
## 26    
## 27    
## 28    
## 29   *
## 30    
## 31    
## 32    
## 33    
## 34    
## 35    
## 36    
## 37   *
## 38    
## 39    
## 40    
## 41    
## 42    
## 43    
## 44    
## 45    
## 46    
## 47    
## 48    
## 49    
## 50    
## 51    
## 52    
## 53    
## 54    
## 55   *
## 56   *

所建立的模型通过回归诊断，影响分析存在异常点，但没有理由怀疑是异常点所以予以保留。

多重共线性

指线性回归模型中的解释变量之间由于存在精确相关关系或高度相关关系而使模型估计失真或难以估计准确。一般来说，由于数据的限制使得模型设计不当，导致设计矩阵中解释变量间存在普遍的相关关系目前常用的多重共线性诊断方法 1.自变量的相关系数矩阵R诊断法：研究变量的两两相关分析，如果自变量间的二元相关系数值很大，则认为存在多重共线性。但无确定的标准判断相关系数的大小与共线性的关系。有时，相关系数值不大，也不能排除多重共线性的可能。

2.方差膨胀因子（the variance inflation factor，VIF)诊断法：方差膨胀因子表达式为：\(VIF_{i}=1/（1-R^{2}_{i})\)。其中Ri为自变量\(x_{i}\)对其余自变量作回归分析的复相关系数。当\(VIF_{i}\)很大时，表明自变量间存在多重共线性。该诊断方法也存在临界值不易确定的问题，在应用时须慎重。

3.容忍值（Tolerance，简记为Tol）法：容忍值实际上是VIF的倒数，即Tol＝1/VIF。其取值在0～1之间，Tol越接近1，说明自变量间的共线性越弱。在应用时一般先预先指定一个Tol值，容忍值小于指定值的变量不能进入方程，从而保证进入方程的变量的相关系数矩阵为非奇异阵，计算结果具有稳定性。但是，有的自变量即使通过了容忍性检验进入方程，仍可导致结果的不稳定。

4.多元决定系数值诊断法：假定多元回归模型p个自变量，其多元决定系数为R2y（X1，X2,…，Xp）。分别构成不含其中某个自变量（\(x_{i}\),i=1,2,…，p）的p个回归模型，并应用最小二乘法准则拟合回归方程，求出它们各自的决定系数\(R^{2}_{i}（i=1,2,…，p）\)。如果其中最大的一个R2k与R2Y很接近，就表明该自变量在模型中对多元决定系数的影响不大，说明该变量对Y总变异的解释能力可由其他自变量代替。它很有可能是其他自变量的线性组合。因此，该自变量进入模型后就有可能引起多重共线性问题。该方法也存在临界值和主观判断问题。

5.条件数与特征分析法：在自变量的观测值构成的设计矩阵X中，求出变量相关系数R的特征值，如果某个特征值很小（如小于0．05 ），或所有特征值的倒数之和为自变量数目的5倍以上，表明自变量间存在多重共线性关系。利用主成分分析，如果X′X的特征值RK小于0．05时，RK所对应的主成分FK可近似为零，表明自变量间存在K个多重共线性关系。

从实际经验的角度,一般若条件数<100,则认为多重共线性的程度很小,若100<=条件数<=1000,则认为存在中等程度的多重共线性,若条件数>1000,则认为存在严重的多重共线性。kappa大于1000，或vif大于10说明存在多重共线性。在R中判断多重共线性的命令为kappa（条件数），vif（方差膨胀因子）

kappa(cor(ex))
## [1] 7.949624
vif(fit) #car包
##        age        sbp log10(wbc) 
##   1.057893   1.012549   1.046398
sqrt(vif(fit)) > 2
##        age        sbp log10(wbc) 
##      FALSE      FALSE      FALSE

一般来说kappa大于1000，或vif大于10说明存在多重共线性,vif开平方是否大于2，若大于2，则存在多重共线性问题。本例中未发现存在多重共线性。

模型比较

AIC(Akaike Information Criterion,赤池信息准则)也可以用来比较模型,它考虑了模型的统计拟合度以及用来拟合的参数数目。AIC值越小的模型要优先选择,它说明模型用较少的参数获得了足够的拟合度。

fit2 <- lm(mtfc~age+sbp+wbc)
AIC(fit,fit2)
##      df      AIC
## fit   5 355.7778
## fit2  5 356.5946

选择AIC值较小的模型lm(mtfc~age+sbp+log10(wbc))。

尽管log10(wbc)回归系数最高，但并不代表log10(wbc)对mtfc的影响最大，因为三个变量的单位不同。对于不同的单位，如果要衡量对因变量大小的影响，需采用标准化回归系数。

lm.beta <- function(MOD) 
{
  b <- summary(MOD)$coef[-1, 1]
  sx <- sapply(MOD$model[-1], sd)
  sy <- sapply(MOD$model[1], sd)
  beta <- b * sx/sy
  return(beta)
}
lm.beta(fit)
##        age        sbp log10(wbc) 
##  0.2375385  0.5440901  0.2404413
detach(ex)

尽管通过模型比较，获得了lm(mtfc~age+sbp+log10(wbc))模型，模型的回归诊断也能通过，但从决定系数来看自变量对因变量的解释程度并不高。

逐步回归

实际中，影响因变量的自变量较多，对如何从自变量中选择若干个，得到最佳的回归方程，在不同的准则下有不同的方法来获得最佳回归方程。对于一个包含n个自变量的的回归问题，全部的回归模型将有\(2^{n}-1\)个。常用的逐步方法有“向前法”，“向后法”，“逐步法”和“最优子集法”。在R中，通过step()函数的direction = c(“both”, “backward”, “forward”)选项分别完成“逐步回归法”、“向后法”和“向前法”。leaps包可以完成全子集回归法,leaps()函数以\(C_{p}\)准则（默认）、校正\(R^{2}\)和\(R^{2}\)来选择全局最优模型。

例 有5个自变量x1～x5和1个因变量，请完成自变量的筛选。

x1<- c(7,1,11,11,7,11,3,1,2,21,1,11,10)
x2<- c(26,29,56,31,52,55,71,31,54,47,40,66,68)
x3<- c(6,15,8,8,6,9,17,22,18,4,23,9,8)
x4<- c(60,52,20,47,33,22,6,44,22,26,34,12,12)
y<- c(78.5,74.3,104.3,87.6,95.9,109.2,102.7,72.5,93.1,115.9,83.8,113.3,109.4)
df <- as.data.frame(cbind(x1,x2,x3,x4,y))

leapmodels <- leaps(x = cbind(x1,x2,x3,x4),y = y)
plot(leapmodels$size, leapmodels$Cp)
abline(0,1)

cbind(leapmodels$size,leapmodels$which, leapmodels$Cp)
##     1 2 3 4           
## 1 2 0 0 0 1 138.730833
## 1 2 0 1 0 0 142.486407
## 1 2 1 0 0 0 202.548769
## 1 2 0 0 1 0 315.154284
## 2 3 1 1 0 0   2.678242
## 2 3 1 0 0 1   5.495851
## 2 3 0 0 1 1  22.373112
## 2 3 0 1 1 0  62.437716
## 2 3 0 1 0 1 138.225920
## 2 3 1 0 1 0 198.094653
## 3 4 1 1 0 1   3.018233
## 3 4 1 1 1 0   3.041280
## 3 4 1 0 1 1   3.496824
## 3 4 0 1 1 1   7.337474
## 4 5 1 1 1 1   5.000000

选择\(C_{p}\)统计量最小的变量集合，本例中\(C_{p}\)统计量最小值所对应的变量集合为x1和x2。结果中1为选中，0为未选中。

subsets <- regsubsets(y~x1+x2+x3+x4,data=df)
summary(subsets)
## Subset selection object
## Call: regsubsets.formula(y ~ x1 + x2 + x3 + x4, data = df)
## 4 Variables  (and intercept)
##    Forced in Forced out
## x1     FALSE      FALSE
## x2     FALSE      FALSE
## x3     FALSE      FALSE
## x4     FALSE      FALSE
## 1 subsets of each size up to 4
## Selection Algorithm: exhaustive
##          x1  x2  x3  x4 
## 1  ( 1 ) " " " " " " "*"
## 2  ( 1 ) "*" "*" " " " "
## 3  ( 1 ) "*" "*" " " "*"
## 4  ( 1 ) "*" "*" "*" "*"
plot(subsets,scale="adjr2") #基于调整R平方，不同子集大小的最佳模型

图的顶部的图形便是最适合的模型，校正R平方值0.98最高，x1和x2两预测变量是最佳模型。也可以通过\(C_{p}\)统计量完成变量的选择。

sbs<- regsubsets(y~x1+x2+x3+x4,data=df)
subsets(sbs,legend=FALSE,statistic="cp",main="cp plot for all subsets regression")
##    Abbreviation
## x1           x1
## x2           x2
## x3           x3
## x4           x4
abline(1,1,lty=2,col="red") #画截距项和斜率均为1的直线

summary(sbs)
## Subset selection object
## Call: regsubsets.formula(y ~ x1 + x2 + x3 + x4, data = df)
## 4 Variables  (and intercept)
##    Forced in Forced out
## x1     FALSE      FALSE
## x2     FALSE      FALSE
## x3     FALSE      FALSE
## x4     FALSE      FALSE
## 1 subsets of each size up to 4
## Selection Algorithm: exhaustive
##          x1  x2  x3  x4 
## 1  ( 1 ) " " " " " " "*"
## 2  ( 1 ) "*" "*" " " " "
## 3  ( 1 ) "*" "*" " " "*"
## 4  ( 1 ) "*" "*" "*" "*"

\(C_{p}\)图越好的模型离截距项和斜率均为1的直线越近，x1-x2、x1-x2-x4和x1-x2-x3-x4均与直线比较接近，这三个模型根据\(C_{p}\)统计量结果类似。

fit <- lm(y~x1+x2+x3+x4,data=df)
fit.step <- step(fit)  #基于AIC
## Start:  AIC=26.94
## y ~ x1 + x2 + x3 + x4
## 
##        Df Sum of Sq    RSS    AIC
## - x3    1    0.1091 47.973 24.974
## - x4    1    0.2470 48.111 25.011
## - x2    1    2.9725 50.836 25.728
## <none>              47.864 26.944
## - x1    1   25.9509 73.815 30.576
## 
## Step:  AIC=24.97
## y ~ x1 + x2 + x4
## 
##        Df Sum of Sq    RSS    AIC
## <none>               47.97 24.974
## - x4    1      9.93  57.90 25.420
## - x2    1     26.79  74.76 28.742
## - x1    1    820.91 868.88 60.629
summary(fit.step)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x4, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0919 -1.8016  0.2562  1.2818  3.8982 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  71.6483    14.1424   5.066 0.000675 ***
## x1            1.4519     0.1170  12.410 5.78e-07 ***
## x2            0.4161     0.1856   2.242 0.051687 .  
## x4           -0.2365     0.1733  -1.365 0.205395    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.309 on 9 degrees of freedom
## Multiple R-squared:  0.9823, Adjusted R-squared:  0.9764 
## F-statistic: 166.8 on 3 and 9 DF,  p-value: 3.323e-08

step()函数通过AIC信息准则，删除了x3变量后AIC值最小24.97，得到y ~ x1 + x2 + x4。回归系数的显著性检验水平有较大提升，但x2和x4的系数检验仍不理想。

drop1(fit.step)
## Single term deletions
## 
## Model:
## y ~ x1 + x2 + x4
##        Df Sum of Sq    RSS    AIC
## <none>               47.97 24.974
## x1      1    820.91 868.88 60.629
## x2      1     26.79  74.76 28.742
## x4      1      9.93  57.90 25.420

去掉x4后，AIC值会上升到25.42，残差的平方和会上升到9.93，是上升最少的。去掉x4后

lm.opt<-lm(y ~ x1+x2, data=df)
summary(lm.opt)
## 
## Call:
## lm(formula = y ~ x1 + x2, data = df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.893 -1.574 -1.302  1.363  4.048 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 52.57735    2.28617   23.00 5.46e-10 ***
## x1           1.46831    0.12130   12.11 2.69e-07 ***
## x2           0.66225    0.04585   14.44 5.03e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.406 on 10 degrees of freedom
## Multiple R-squared:  0.9787, Adjusted R-squared:  0.9744 
## F-statistic: 229.5 on 2 and 10 DF,  p-value: 4.407e-09
gvlma(lm.opt)
## 
## Call:
## lm(formula = y ~ x1 + x2, data = df)
## 
## Coefficients:
## (Intercept)           x1           x2  
##     52.5773       1.4683       0.6623  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = lm.opt) 
## 
##                        Value p-value                Decision
## Global Stat        1.5244265  0.8223 Assumptions acceptable.
## Skewness           0.5484212  0.4590 Assumptions acceptable.
## Kurtosis           0.5559113  0.4559 Assumptions acceptable.
## Link Function      0.0004057  0.9839 Assumptions acceptable.
## Heteroscedasticity 0.4196883  0.5171 Assumptions acceptable.
influence.measures(lm.opt)
## Influence measures of
##   lm(formula = y ~ x1 + x2, data = df) :
## 
##     dfb.1_  dfb.x1  dfb.x2  dffit cov.r cook.d    hat inf
## 1  -0.3880 -0.0625  0.3558 -0.428 1.535 0.0639 0.2512    
## 2   0.2597 -0.1375 -0.1649  0.291 1.719 0.0305 0.2619    
## 3   0.0563 -0.1000 -0.0754 -0.239 1.356 0.0202 0.1189    
## 4  -0.2993 -0.2256  0.3277 -0.439 1.491 0.0668 0.2422    
## 5  -0.0127  0.0244 -0.0480 -0.176 1.339 0.0111 0.0836    
## 6  -0.1410  0.3248  0.1961  0.742 0.488 0.1387 0.1151    
## 7   0.2147  0.2671 -0.4008 -0.495 1.867 0.0867 0.3618    
## 8  -0.4816  0.2850  0.2866 -0.558 1.327 0.1039 0.2412    
## 9   0.0288 -0.2733  0.1584  0.385 1.344 0.0510 0.1792    
## 10 -0.0239  0.8515 -0.2205  0.919 2.442 0.2903 0.5500   *
## 11  0.4911 -0.5201 -0.1450  0.769 0.781 0.1695 0.1840    
## 12 -0.1058  0.0430  0.1278  0.189 1.627 0.0131 0.1967    
## 13  0.4379 -0.0667 -0.5607 -0.744 0.949 0.1672 0.2142
sqrt(vif(lm.opt)) > 2
##    x1    x2 
## FALSE FALSE

检验结果显著，回归诊断和多重共线性均通过，影响分析中只发现有一个异常点。因为全子集回归考虑了更多模型，全子集回归要优于逐步回归，但当自变量个数较多时，全子集回归较慢。一般来所，变量的自动筛选应建立在背景知识理解基础上进行，防止出现拟合效果好，但没有实际意义的模型。

交叉验证

就是按一定比例将原始数据按照拆分成训练集和测试集，现在训练集上获取回归方程，然后在测试集上做预测。由于测试集不涉及模型参数的选择，该样本可获得比新数据获得更为精确的估计。\(k\)重交叉验证中，样本被分为\(k\)个子样本，轮流将\(k-1\)个子样本作为训练集，另外1个样本作为测试集。通过获得\(k\)个预测方程，记录\(k\)个测试集的预测结果，然后求其平均值。

shrinkage<-function(fit,k=5){  
  require(bootstrap)  
  theta.fit<-function(x,y){lsfit(x,y)}  
  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-Validated R-square=",r2cv,"\n")  
  cat("Change=",r2-r2cv,"\n")  
}  

fit <- lm(y~x1+x2,data=df)
shrinkage(fit)
## Original R-square= 0.9786784 
## 5 Fold Cross-Validated R-square= 0.9664616 
## Change= 0.01221677

获得原始R平方为0.9786,交叉验证后的R平方为0.9962（基于BootStrap方法，每次运行结果会有不同）。R平方减少得越少，预测越精准。

相对重要性

评价自变量相对重要性，最简单的方法为比较标准化的回归系数，它表示当其他自变量不变时，该自变量变化1个单位引起的因变量的变化。前面通过lm.beta()函数获得了标准化的回归系数。基于相对权重的重要性测量，是对所有可能自模型添加一个自变量引起的R平方平均增加量的一个近似值，比标准化回归系数更为直观。

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)<-"Weight"  
  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)  
}  
fit <- lm(y~x1+x2,data=df)
relweights(fit,col="lightgrey")  

##      Weight
## x1 43.23985
## x2 56.76015

可以看到x2解释了56.7%的R平方，x1解释了43.2的平方，x2相比x1更为重要。

分位数回归

传统的线性回归模型描述了因变量的条件均值分布受自变量的影响过程。最小二乘法是估计回归系数的最常用的方法。如果模型的随机误差项来自均值为零、方差相同的分布，那么模型回归系数的最小二乘估计为最佳线性无偏估计（BLUE）；如果随机误差项是正态分布，那么模型回归系数的最小二乘估计与极大似然估计一致，均为最小方差无偏估计（MVUL）。分位数回归(Quantile Regression)利用解释变量的多个分位数（例如四分位、十分位、百分位等）来得到被解释变量的条件分布的相应的分位数方程。与传统的OLS只得到均值方程相比，它可以更详细地描述变量的统计分布。 在数据出现尖峰或厚尾的分布、存在显著的异方差等情况，传统的线性回归模型的假设常常不被满足，最小二乘法估计将不再具有上述优良性且稳健性非常差。最小二乘回归假定自变量只能影响因变量的条件分布的位置，但不能影响其分布的刻度或形状的任何其他方面。分位数回归依据因变量的条件分位数对自变量进行回归，这样得到了所有分位数下的回归模型。因此分位数回归相比普通最小二乘回归只能描述自变量对于因变量局部变化的影响而言，更能精确地描述自变量对于因变量的变化范围以及条件分布形状的影响。分位数回归能够捕捉分布的尾部特征，当自变量对不同部分的因变量的分布产生不同的影响时．例如出现左偏或右偏的情况时。它能更加全面的刻画分布的特征，从而得到全面的分析，而且其分位数回归系数估计比OLS回归系数估计更稳健。

例 quantreg包中自带数据集engel描述了食物支出与家庭收入之间关系，其数据格式如下，请完成分位数回归。

data(engel,package = "quantreg")
pander(head(engel))
income foodexp
420.2 255.8
541.4 311
901.2 485.7
639.1 403
750.9 495.6
945.8 633.8 
#进行分位数回归
fit = rq(foodexp ~ income, tau = c(0.1,0.25,0.5,0.75,0.9), 
         data = engel,method = "br")  
summary(fit)
## 
## Call: rq(formula = foodexp ~ income, tau = c(0.1, 0.25, 0.5, 0.75, 
##     0.9), data = engel, method = "br")
## 
## tau: [1] 0.1
## 
## Coefficients:
##             coefficients lower bd  upper bd 
## (Intercept) 110.14157     79.88753 146.18875
## income        0.40177      0.34210   0.45079
## 
## Call: rq(formula = foodexp ~ income, tau = c(0.1, 0.25, 0.5, 0.75, 
##     0.9), data = engel, method = "br")
## 
## tau: [1] 0.25
## 
## Coefficients:
##             coefficients lower bd  upper bd 
## (Intercept)  95.48354     73.78608 120.09847
## income        0.47410      0.42033   0.49433
## 
## Call: rq(formula = foodexp ~ income, tau = c(0.1, 0.25, 0.5, 0.75, 
##     0.9), data = engel, method = "br")
## 
## tau: [1] 0.5
## 
## Coefficients:
##             coefficients lower bd  upper bd 
## (Intercept)  81.48225     53.25915 114.01156
## income        0.56018      0.48702   0.60199
## 
## Call: rq(formula = foodexp ~ income, tau = c(0.1, 0.25, 0.5, 0.75, 
##     0.9), data = engel, method = "br")
## 
## tau: [1] 0.75
## 
## Coefficients:
##             coefficients lower bd  upper bd 
## (Intercept)  62.39659     32.74488 107.31362
## income        0.64401      0.58016   0.69041
## 
## Call: rq(formula = foodexp ~ income, tau = c(0.1, 0.25, 0.5, 0.75, 
##     0.9), data = engel, method = "br")
## 
## tau: [1] 0.9
## 
## Coefficients:
##             coefficients lower bd  upper bd 
## (Intercept)  67.35087     37.11802 103.17399
## income        0.68630      0.64937   0.74223
# 通过设置参数se，可以得到系数的假设检验
summary(fit, se = "nid")  
## 
## Call: rq(formula = foodexp ~ income, tau = c(0.1, 0.25, 0.5, 0.75, 
##     0.9), data = engel, method = "br")
## 
## tau: [1] 0.1
## 
## Coefficients:
##             Value     Std. Error t value   Pr(>|t|) 
## (Intercept) 110.14157  29.39768    3.74661   0.00023
## income        0.40177   0.04024    9.98420   0.00000
## 
## Call: rq(formula = foodexp ~ income, tau = c(0.1, 0.25, 0.5, 0.75, 
##     0.9), data = engel, method = "br")
## 
## tau: [1] 0.25
## 
## Coefficients:
##             Value    Std. Error t value  Pr(>|t|)
## (Intercept) 95.48354 21.39237    4.46344  0.00001
## income       0.47410  0.02906   16.31729  0.00000
## 
## Call: rq(formula = foodexp ~ income, tau = c(0.1, 0.25, 0.5, 0.75, 
##     0.9), data = engel, method = "br")
## 
## tau: [1] 0.5
## 
## Coefficients:
##             Value    Std. Error t value  Pr(>|t|)
## (Intercept) 81.48225 19.25066    4.23270  0.00003
## income       0.56018  0.02828   19.81032  0.00000
## 
## Call: rq(formula = foodexp ~ income, tau = c(0.1, 0.25, 0.5, 0.75, 
##     0.9), data = engel, method = "br")
## 
## tau: [1] 0.75
## 
## Coefficients:
##             Value    Std. Error t value  Pr(>|t|)
## (Intercept) 62.39659 16.30538    3.82675  0.00017
## income       0.64401  0.02324   27.71244  0.00000
## 
## Call: rq(formula = foodexp ~ income, tau = c(0.1, 0.25, 0.5, 0.75, 
##     0.9), data = engel, method = "br")
## 
## tau: [1] 0.9
## 
## Coefficients:
##             Value    Std. Error t value  Pr(>|t|)
## (Intercept) 67.35087 22.39538    3.00736  0.00292
## income       0.68630  0.02849   24.08853  0.00000
plot(fit)

tau表示计算多个分位点的分位数回归结果，如tau = c(0.25,0.5,0.75)是同时计算25%、50%、75%分位数下的回归结果。method：进行拟合的方法，取值包括：默认值“br”，表示 Barrodale & Roberts算法的修改版；“fn”，针对大数据可以采用的Frisch–Newton内点算法；“pfn”，针对特别大数据，使用经过预处理的Frisch–Newton逼近方法；“fnc”，针对被拟合系数特殊的线性不等式约束情况； “lasso”和“scad”，基于特定惩罚函数的平滑算法进行拟合。 se = “rank”: 按照Koenker(1994)的排秩方法计算得到的置信区间，默认残差为独立同分布。注意的是，上下限是不对称的。se=”iid”: 假设残差为独立同分布，用KB（1978）的方法计算得到近似的协方差矩阵。se = “nid”: 表示按照Huber方法逼近得到的估计量。se=”ker”:采用Powell(1990)的核估计方法。se=”boot”:采用bootstrap方法自助抽样的方法估计系数的误差标准差。

穷人和富人的消费比较

data(engel,package = "quantreg")
attach(engel)
z=rq(foodexp~income,tau=-1,data = engel) #tau不再[0,1]时，表示按最细分位点划分
x.poor=quantile(income,0.1) #10%分位点的收入，穷人
x.rich=quantile(income,0.9) #90%分位点的收入，富人
ps=z$sol[1,] #每个分位点的tau值
qs.poor=c(c(1,x.poor)%*%z$sol[4:5,]) #穷人的消费估计值
qs.rich=c(c(1,x.rich)%*%z$sol[4:5,]) #富人的消费估计值
par(mfrow=c(1,2))
plot(c(ps,ps),c(qs.poor,qs.rich),type="n",     # type=”n”表示初始化图形区域，但不画图
     xlab=expression(tau), ylab="quantile")
plot(stepfun(ps,c(qs.poor[1],qs.poor)), do.points=F,
     add=T)
plot(stepfun(ps,c(qs.poor[1],qs.rich)), do.points=F,
     add=T, col.hor="gray", col.vert="gray")

ps.wts = ( c(0,diff(ps)) + c(diff(ps),0) )/2
ap = akj(qs.poor, z=qs.poor, p=ps.wts)
ar = akj(qs.rich, z=qs.rich, p=ps.wts)
plot(c(qs.poor,qs.rich), c(ap$dens, ar$dens),
     type="n", xlab="Food Expenditure", ylab="Density")
lines(qs.rich,ar$dens,col="gray")
lines(qs.poor,ap$dens,col="black")
legend("topright", c("poor","rich"), lty=c(1,1),
       col=c("black","gray"))

par(mfrow=c(1,1))

穷人和富人的食品消费支出有明显的不同，穷人在不同分位点食品消费支出差别不大，富人不同分位点食品消费支出差别较大。右图表示，穷人消费支出集中于400左右，富人消费支出集中于800～1200。

模型比较

# 比较不同分位点下，收入对食品支出的影响机制是否相同
fit1 = rq(foodexp ~ income, tau = 0.25, data = engel)
fit2 = rq(foodexp ~ income, tau = 0.5, data = engel)
fit3 = rq(foodexp ~ income, tau = 0.75, data = engel)
anova(fit1,fit2,fit3)
## Quantile Regression Analysis of Deviance Table
## 
## Model: foodexp ~ income
## Joint Test of Equality of Slopes: tau in {  0.25 0.5 0.75  }
## 
##   Df Resid Df F value    Pr(>F)    
## 1  2      703  15.557 2.449e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

P值远小于0.05，故不同分位点下收入对食品支出的影响机制不同。 不同分位点拟合曲线比较

plot(income,foodexp,cex=0.25,type = "n",
     xlab = "Household Income",ylab = "Food Expenditure")
points(income,foodexp,cex=0.5,col="blue")
abline(rq(foodexp~income,tau = 0.5),col="blue")
abline(lm(foodexp~income),lty=2,col="red")
taus=c(0.1,0.25,0.75,0.9)
for(i in 1:length(taus)){
  abline(rq(foodexp~income,tau=taus[i]),col="gray")
}

detach(engel)

残差形态检验

#位值漂移模型：不同分位点估计结果之间的斜率相同或相近，截距不同
KhmaladzeTest(foodexp ~ income, 
              data = engel, taus = seq(.05,.95,by = .01),nullH = "location")
## taus: 0.05  0.06  0.07  0.08  0.09  0.1  0.11  0.12  0.13  0.14  0.15  0.16  0.17  0.18  0.19  0.2  0.21  0.22  0.23  0.24  0.25  0.26  0.27  0.28  0.29  0.3  0.31  0.32  0.33  0.34  0.35  0.36  0.37  0.38  0.39  0.4  0.41  0.42  0.43  0.44  0.45  0.46  0.47  0.48  0.49  0.5  0.51  0.52  0.53  0.54  0.55  0.56  0.57  0.58  0.59  0.6  0.61  0.62  0.63  0.64  0.65  0.66  0.67  0.68  0.69  0.7  0.71  0.72  0.73  0.74  0.75  0.76  0.77  0.78  0.79  0.8  0.81  0.82  0.83  0.84  0.85  0.86  0.87  0.88  0.89  0.9  0.91  0.92  0.93  0.94  0.95
## $nullH
## [1] "location"
## 
## $Tn
## [1] 1.515467
## 
## $THn
## [1] 1.515467
## 
## attr(,"class")
## [1] "KhmaladzeTest"
KhmaladzeTest(foodexp ~ income, 
              data = engel, taus = seq(.05,.95,by = .01),nullH = "location",se="ker")
## taus: 0.05  0.06  0.07  0.08  0.09  0.1  0.11  0.12  0.13  0.14  0.15  0.16  0.17  0.18  0.19  0.2  0.21  0.22  0.23  0.24  0.25  0.26  0.27  0.28  0.29  0.3  0.31  0.32  0.33  0.34  0.35  0.36  0.37  0.38  0.39  0.4  0.41  0.42  0.43  0.44  0.45  0.46  0.47  0.48  0.49  0.5  0.51  0.52  0.53  0.54  0.55  0.56  0.57  0.58  0.59  0.6  0.61  0.62  0.63  0.64  0.65  0.66  0.67  0.68  0.69  0.7  0.71  0.72  0.73  0.74  0.75  0.76  0.77  0.78  0.79  0.8  0.81  0.82  0.83  0.84  0.85  0.86  0.87  0.88  0.89  0.9  0.91  0.92  0.93  0.94  0.95
## $nullH
## [1] "location"
## 
## $Tn
## [1] 1.335624
## 
## $THn
## [1] 1.335624
## 
## attr(,"class")
## [1] "KhmaladzeTest"
#位置-尺度漂移模型：不同分位点估计结果斜率和截距都不同
KhmaladzeTest(foodexp ~ income, 
              data = engel, taus = seq(.05,.95,by = .01),nullH = "location-scale")
## taus: 0.05  0.06  0.07  0.08  0.09  0.1  0.11  0.12  0.13  0.14  0.15  0.16  0.17  0.18  0.19  0.2  0.21  0.22  0.23  0.24  0.25  0.26  0.27  0.28  0.29  0.3  0.31  0.32  0.33  0.34  0.35  0.36  0.37  0.38  0.39  0.4  0.41  0.42  0.43  0.44  0.45  0.46  0.47  0.48  0.49  0.5  0.51  0.52  0.53  0.54  0.55  0.56  0.57  0.58  0.59  0.6  0.61  0.62  0.63  0.64  0.65  0.66  0.67  0.68  0.69  0.7  0.71  0.72  0.73  0.74  0.75  0.76  0.77  0.78  0.79  0.8  0.81  0.82  0.83  0.84  0.85  0.86  0.87  0.88  0.89  0.9  0.91  0.92  0.93  0.94  0.95
## $nullH
## [1] "location-scale"
## 
## $Tn
## [1] 0.7097263
## 
## $THn
## [1] 0.7097263
## 
## attr(,"class")
## [1] "KhmaladzeTest"
KhmaladzeTest(foodexp ~ income, 
              data = engel, taus = seq(.05,.95,by = .01),nullH = "location-scale",se="ker")
## taus: 0.05  0.06  0.07  0.08  0.09  0.1  0.11  0.12  0.13  0.14  0.15  0.16  0.17  0.18  0.19  0.2  0.21  0.22  0.23  0.24  0.25  0.26  0.27  0.28  0.29  0.3  0.31  0.32  0.33  0.34  0.35  0.36  0.37  0.38  0.39  0.4  0.41  0.42  0.43  0.44  0.45  0.46  0.47  0.48  0.49  0.5  0.51  0.52  0.53  0.54  0.55  0.56  0.57  0.58  0.59  0.6  0.61  0.62  0.63  0.64  0.65  0.66  0.67  0.68  0.69  0.7  0.71  0.72  0.73  0.74  0.75  0.76  0.77  0.78  0.79  0.8  0.81  0.82  0.83  0.84  0.85  0.86  0.87  0.88  0.89  0.9  0.91  0.92  0.93  0.94  0.95
## $nullH
## [1] "location-scale"
## 
## $Tn
## [1] 0.6224596
## 
## $THn
## [1] 0.6224596
## 
## attr(,"class")
## [1] "KhmaladzeTest"

Tn表示模型整体检验，THn表示每个自变量的检验。位值漂移模型的Tn值比位置-尺度漂移模型的Tn值大，拒绝位值漂移模型的概论较大，位值-尺度漂移模型更加合适。

分位数回归的分解

分位数分解法对各个影响因素进行分解分析

# MM2005分位数分解的函数
MM2005 = function(formu,taus, data, group, pic=F){
  # furmu 为方程，如foodexp~income
  # taus 为不同的分位数
  # data 总的数据集
  # group 分组指标，是一个向量，用于按行区分data
  # pic 是否画图，如果分位数比较多，建议不画图
  engel1 = data[group==1,]
  engel2 = data[group==2,]
  # 开始进行分解
  fita = summary( rq(formu, tau = taus, data = engel1 ) )
  fitb = summary( rq(formu, tau = taus, data = engel2 ) )
  tab = matrix(0,length(taus),4)
  colnames(tab) = c("分位数","总差异","回报影响","变量影响")
  rownames(tab) = rep("",dim(tab)[1])
  for( i in 1:length(taus) ){
    ya = cbind(1,engel1[,names(engel1)!=formu[[2]]] ) %*% fita[[i]]$coef[,1]
    yb = cbind(1,engel2[,names(engel2)!=formu[[2]]] ) %*% fitb[[i]]$coef[,1]
    # 这里以group==1为基准模型，用group==2的数据计算反常规模型拟合值
    ystar = cbind(1,engel2[,names(engel2)!=formu[[2]]] ) %*% fita[[i]]$coef[,1]
    ya = mean(ya)
    yb = mean(yb)
    ystar = mean(ystar)
    tab[i,1] = fita[[i]]$tau
    tab[i,2] = yb - ya
    tab[i,3] = yb - ystar # 回报影响，数据相同，模型不同：模型机制的不同所产生的差异
    tab[i,4] = ystar - ya # 变量影响，数据不同，模型相同：样本点不同产生的差异
  }
  # 画图
  if( pic ){
    attach(engel)
    windows(5,5)
    plot(income, foodexp, cex=0.5, type="n",main="两组分位数回归结果比较")
    points(engel1, cex=0.5, col=2)
    points(engel2, cex=0.5, col=3)
    for( i in 1:length(taus) ){
      abline( fita[[i]], col=2 )
      abline( fitb[[i]], col=3 )
    }
    detach(engel)
  }
  # 输出结果
  tab
}

data(engel,package = "quantreg")
group = c(rep(1,100),rep(2,135))  # 取前100个为第一组，后135个第二组
taus = c(0.05,0.25,0.5,0.75,0.95)  # 需要考察的不同分位点
MM2005(foodexp~income, taus, data = engel, group=group, pic=F)
##  分位数     总差异  回报影响 变量影响
##    0.05 -30.452061 -72.35939 41.90733
##    0.25  -2.017317 -46.20125 44.18394
##    0.50  30.941212 -23.24042 54.18163
##    0.75  43.729025 -15.76283 59.49185
##    0.95  52.778985 -11.29932 64.07830