直线回归可以用以下直线回归方程（linear regression equation）来表示，如下所示：
公式（一）称为经验回归方程或样本回归方程，其中b表示这条方程的斜率，a表示这个方程在y轴上的截矩，有的教材上把a与b写作是$\hat{a}$和$\hat{b}$，这两个字母上面的符号读作hat，称为“帽子”，它们表示利用样本的数据估计得得来的截矩和斜率。我们可以通过一系列的计算求出这个方程的a和b，这个方程是对两变量总体间线性关系的一个估计，根据散点图可以假设，对于X的各个取值，相应Y的总体均数$\mu_{Y|X}$在一条直线上，如下所示：
总体均数表示为：
除了图中所示两变量呈直线关系外，一般还假定每个X对应Y的总体为正态分布，各个正态分布的总体方差相等且各次观测相互独立，这样公式（一）中的$\hat{Y}$实际上是x所对应Y的总体均数$\mu_{Y|X}$的一个样本估计值，称为回归方程的预测值（predicted value），而a、 b分别为$\alpha$和$\beta$的样本估计，其中a称为常数项，b称为回归系数（coefficient of regression），b是直线的斜率（slope），其统计意义是，当X变化一个单位时，Y的平均改变的估计值，b>0时，直线从左下方走到右上方，Y随X的增大而增大，当b<0时，直线从左上方走向右下方，Y随X的增大而减小；当b=0时，直线与X轴平行，Y与X没有关系。
直线回归方程的计算思路
如果能够从样本数据中求得a和b的数值，那么回归方程即可唯一确定，从散点图上来看，求解a和b实际上就是怎么找到一条最能代表数据点分布趋势的直线，将实测值Y与假定回归线上估计值$\hat{Y}$（这个Y的上面有一个帽子结构的符号表示这是Y值的估计值）的纵向距离$Y-\hat{Y}$称为残差（residual），残差通常用$e_{i}$表示，如果把残差进行了均一化，即进行了z转换，那它就是标准化残差，其中对于要找到的这条最佳的直线，它的残差平方和（sum of squared residuals, sum of squared errors或residual sum of squares，简称为SSE或RSS等）是最小的，这就是最小二乘法（Least sum of squares，LS）的思想（因为古中国的数学家在描述平方时，使用的术语是“二乘”，因此这里叫最小二乘法），如下所示：
在一定假设条件下，如此得到的回归系数最为理解，按照这一原则，数学上可以很容易得到a和b的计算公式，如下所示：
$b=\frac{l_{XY}}{l_{XX}}=\frac{\sum(X-\bar{X})(Y-\bar{Y}))}{\sum(X-\bar{X})^2}$ 公式（三）
$a=\bar{Y}-b \bar{X}$ 公式（四）
其中，$l_{XY}$是X与Y的均离均差平交叉乘积和，简称离均差积和，公式为：
$l_{XY}=\sum(X-\bar{X})(Y-\bar{Y})=\sum XY-\frac{\sum(X)(\sum(Y)}{n}$ 公式（五）
除了用公式（一）来表示两变量线性回归关系，还可以在散点图上绘制出样本回归直线作为一种直
观的统计描述补充形式，此直线必然通过点$(\bar{X},\bar{Y})$，且与纵坐标轴相关于截矩a，如果散点图没有从坐标系原点开始，可以在自变量实测范围内远端取易于读取的X值代入回归方程得到一个点的坐标，连接此点与点$(\bar{X},\bar{Y})$也能绘制出回归直线。
现在再看一下原始数据：
| 编号 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
|---|---|---|---|---|---|---|---|---|
| 年龄X | 13 | 11 | 9 | 6 | 8 | 10 | 12 | 7 |
| 尿肌酐含量Y | 3.54 | 3.01 | 3.09 | 2.48 | 2.56 | 3.36 | 3.18 | 2.65 |
由原始数据与散点图可知，两变量之间呈直线趋势，因此可以计算出相应的参数：
线性回归方程的计算过程
第一步：计算出X、Y的均值$\bar{X}$、$\bar{Y}$，离均差平方和$l_{XX}$、$l_{YY}$与离均差积和$l_{XY}$，如下所示：
$\bar{X}=\frac{\sum X}{n}=\frac{76}{8}=9.5,\bar{Y}=\frac{\sum X}{n}=\frac{23.87}{8}=9.5$
$l_{XX}=\sum X^2-\frac{(\sum X)^2}{n}=764-\frac{(76)^2}{8}=42$
$l_{YY}=\sum X^2-\frac{(\sum Y)^2}{n}=72.2683-\frac{(23.87)^2}{8}=1.0462$
$l_{YY}=\sum XY-\frac{(\sum X)(\sum Y)}{n}=232.61-\frac{(76)(23.87)}{8}=5.8450$
第二步：计算出回归系数b和截矩a
根据公式（三），b=5.8450/42=0.1392
根据公式（四），a=2.9838-(0.1392)(9.5)=1.667
列出直线回归方程，如下所示：
$\hat{Y}=0.1392X + 1.667$
R计算结果如下所示：
|
|
直线回归中的统计推断
回归方程的假设检验
建立样本直线回归方程，只是完成了统计分析中两个变量关系的描述，研究者还必须要回答它所来自总体的直线回归关系是否确实存在，即是否对总体有$\beta \neq 0$（$\beta$就是自变量与因变量的相关系数），如下图所示，无论X如何取值，$\mu_{Y|X}$总在一条水平线上，即$\beta=0$，总体直线回归方程并不成立，也就是说Y与X无直接线性关系，此时$\mu_{Y|X}=\mu_{Y}$，然而在一次随机抽样中，如果所得的样本为实心圆点所示，则会得到一个并不等于0的样本回归系数b，b与0相关多大可以认为具有统计学意义？这就需要方差分析或与其等价的t检验来说明问题。
方差分析
先看一张图，如下所示：
在上图中，任意一点P的纵坐标被回归直线$\hat{Y}$与无数$\bar{Y}$截成三个线段，其中：$Y-\bar{Y}=(\hat{Y}-\bar{Y}=(Y-\hat{Y}))$。由于P点是散点图中任取的一点，将全部数据点都这样处理，并将等式两端平方后再求和，可以证明：
$\sum(\hat{Y}-\bar{Y})(Y-\hat{Y})=0$
那么有以下的公式：
$\sum(Y-\hat{Y})^2=\sum(\hat{Y}-\bar{Y})^2+\sum(\bar{Y}-\hat{Y})^2$
上面的公式用符号表示就是：$SS_{总}=SS_{回}+SS_{残}$（公式六)
其中，$SS_{总}$即$\sum(Y-\bar{Y})^2$，这是Y的离均差平方和，表示未考虑Y与X的回归关系时，Y的总变异。
$SS_{回}$即$\sum(\hat{Y}-\bar{Y})^2$，这是回归平方和，由于特定样本的均数$\bar{Y}$是固定的，因此这部分变异由$\hat{Y}_{i}$的大小不同引起，当X被引物回归以后，正是由于$X_{i}$的不同导致了$\hat{Y}=a+bX_{i}$的不同，所以，$SS_{回}$就反映了在Y的总变异中可以用Y与X的直线关系解释的那部分变异，b离0越远，Y受X的影响越大，$SS_{回}$就越大，说明回归效果就越好。
$SS_{残}$即$\sum({Y-\hat{Y}})^2$，为残差平方和。它反映了除了X对Y的线性影响之外的一切因素对Y的变异的作用，也就是在总平方和中无法用X解释的部分，表示考虑回归之后Y真正的随机误差，在散点图中，各实测点离回归直线越近，$SS_{残}$也就越小，说明直线回归的估计误差越小，回归的作用就越明显。
上述的三个平方和，各有其相应的自由度$\nu$，并有如下的关系：
$\nu_{总}=\nu_{回}+\nu_{残}, \nu_{总}=n-1, \nu_{回}=1, \nu_{残}=n-2$
从以上离均差平方和及其自由度的分解可知，不考虑回归时，随机误差是Y的总变异$SS_{总}$，而考虑回归以后，由于回归的贡献，使原来的随机误差减小为$SS_{残}$，如果两变量之间总体回归关系确实存在，回归的贡献就要大于随机误差，大到何种程度时可以认为具有统计学意义，可以计算如下F统计量（其实就是方差分析了），如下所示：
$F=\frac{SS_{回}/\nu_{回}}{SS_{残}/\nu_{残}}=\frac{MS_{回}}{MS_{残}}, \nu_{回}=1, \nu_{残}=n-2$（公式八）
在这个公式中，$MS_{回}$、$MS_{残}$分别称为回归均方与残差均方，统计量F服从自由度为$\nu_{回}$，$\nu_{残}$的F分布，求F值后，查F界值表，得到p值，按所取检验水准作为推断结论。实际计算时，可以先将$X_{i}$依次代入回归方程，求得${\hat{Y}}$，再求得$SS_{残}$，$SS_{回}$也可以利用下面的公式直接求出$SS_{回}$，再得到$SS_{残}$，如下所示：
$SS_{回}=bl_{XY}=l^2_{XY}/l_{XX}=b^2l_{XX}$（公式九）
t检验
对于$\beta=0$这一假设是否成立还可以进行如下t检验，如下所示：
$t=\frac{b-0}{S_{b}},\nu=n-2$（公式十）
$S_{b}=\frac{S_{Y \cdot X}}{\sqrt{L_{XX}}}$（公式十一）
$S_{X \cdot Y}=\sqrt{\frac{SS_{残}}{n-2}}$（公式十二）
其中，$SS_{Y\cdot X}$是回归的剩余标准差（standard deviation of residuals），$S_{b}$是样本回归系数标准误，这些计算公式表明，扩大自变量的取值范围可以减小$S_{b}$，使得回归系数的估计更稳定。
第1案例：简单线性回归
现在我们看一下前面的数据得到的直线回归方程是否成立，前面得到的线性回归方程如下所示：
$\hat{Y}=0.1392X + 1.667$
方差分析
$H_{0}:\beta=0$，即尿肌酐含量与年龄之间无直线关系。
$H_{1}:\beta \neq 0$，即尿肌酐含量与年龄之间有直线关系。
$\alpha=0.05$
根据公式（九），可知$SS_{回}=l^2_{XY}/l_{XX}=5.845^2/42=0.8134$
根据公式（六），可知$SS_{残}=SS_{总}-SS_{回}=1.0462-0.8134=0.2328$
下面是方差分析表，如下所示：
| 变异来源 | 自由度 | SS | MS | F | P |
|---|---|---|---|---|---|
| 变异来源 | 7 | 10462 | |||
| 回归 | 1 | 0.8134 | 0.8134 | 20.963 | <0.01 |
| 残差 | 6 | 0.2328 | 0.0388 |
$\nu_{1}=1$、$\nu_{2}=6$，查F界值表，得P<0.01，按$\alpha=0.05$的检验水准，拒绝$H_{0}$，接受$H_{1}$，可以认为尿肌酐含量与年龄宰有直线关系。
在Ｒ对这个案例进行方差检验，如下所示：
|
|
t检验
t检验的$H_{0}$、${H_{1}}$与$\alpha$同上，在这个案例中，n=8，$SS_{残}=0.2328$、$L_{XX}=42$，b=0.1392，根据公式（十），公式（十一），公式（十二），可得：
$S_{Y\cdot X}=\sqrt{\frac{0.2328}{8-2}}=0.1970$, S_{b}=\frac{0.1970}{\sqrt{42}}=0.0304$$
$t=\frac{0.1392}{0.0304}=4.579$
由于$\nu=6$，查t界值表可知，0.002<P<0.005，按照$\alpha=0.05$的检验水准，拒绝$H_{0}$，接受$H_{1}$，可以认为尿肌酐含量与年龄间有直线关系。注意，在本例中，$\sqrt{F}=\sqrt{20.97}=4.579=t$，实际上直线回归中对回归系数的t检验与F检验等价，类似于两样本均数比较可以作为t检验，也可以做方差分析。
这里要说明一下，两个变量关系的密切程度或数量上的影响大小的统计量是相关系数或者是另归系数的绝对值r，而不是假设检验的p值，p值越小，只能说越有理由认为变量间的直线关系存在，而不能说关系越密切或越“显著”。另外，直线回归用于预测时，其适用范围一般不应超出样本中自变量的取值范围，此时求得的预测值称内插（interpolation），而超过自变量取值范围所得预测值称为外延（extrapolation），若无充分理由说明现有自变量范围以外的两变量间仍然是直线关系，应尽量避免不合理的外延，举个例子，我们在使用BCA测蛋白浓度时，标准品的浓度范围是0-0.5ug/uL，我们构建的线性回归方程是蛋白浓度=a X OD562 + b，如果我们的样本OD562的吸光度超过了标准品最高浓度的OD562的值（例如标准品0.5ug/uL的OD562值是0.4，我的蛋白样本的OD562的值是0.45），那么此时，我们使用求出来的线性回归方程来计算蛋白样本中的蛋白浓度是不太准的。
现在用R来计算一下，如下所示：
|
|
结果解释：
其中0.4.579是对系数（斜率）是否为零t的检验，t=4.579,p<0.00377，表示斜率不为零，它的零假设为$H_{0}:\beta=0$，而备选假设为$H_{1}:\beta \neq 0$。
Residual standard error: 0.197 on 6 degrees of freedom，这个是残差，自由度为6
Multiple R-squared: 0.7775, Adjusted R-squared: 0.7404，R平方，叫测定系数或可决系数，越接近于1表示两者的线性关系越好。
F-statistic: 20.97 on 1 and 6 DF, p-value: 0.003774，F检验，p值小于0.01，在只一个自变量的情况下，F检验的p值与t检验的p值在数值上相等，并且F的值等于t值的平方。
总体回归系数$\beta$的可信区间
利用上述对回归系数的t检验，我们可以得到$\beta$的$1-\alpha$的置信区间为
$b\pm t_{\alpha/2,\nu}S_{b}$ （公式十三）
根据前面计算的结果，b=0.1392，$S_{b}=0.0304$，按自由度为6，查t界值表，得到t值为2.447，那么根据公式（十三），我们可以计算出$\beta$的95%置信区间，如下所示：
$(0.1392-2.447 \times0.0304,0.1392+2.447 \times0.0304)$
即$(0.0648,0.2136)$，需要注意的是，这个区间不包括0，可按$\alpha=0.05$水准同样得到总体回归系数不为0的结论，即用区间估计回答相同$\alpha$时的假设检验问题。
利用回归方程进行估计和预测
总体均数$\mu_{Y|X}$的可信区间
给定X的数值$X_{0}$，由样本回归方程计算出的$\hat{Y}$只是相应总体均数$\mu_{Y|X}$的一个点估计，$\hat{Y}_{0}$会因样本而异，存在抽样误差，反映其抽样误差的标准误按以下公式计算：
$S_{\hat{Y}_{0}}=S_{Y\cdot X}\sqrt{\frac{1}{n}-\frac{(X_{0}-\bar{X})^2}{\sum(X-\bar{X})^2}}$（公式十四）
$\hat{Y}$的标准误与回归的剩余标准差$S_{Y\cdot X}$成正比，当$X_{0}=\bar{X}$时，达到其最小值$S_{Y\cdot X}/\sqrt{n}$，$X_{0}$离$\bar{X}$越远，其标准误$S_{\hat{Y}_{0}}$越大，给定$X=X_{0}$时，总体均数$\nu_{Y|X_{0}}$的$1-\alpha$的置信区间为
$\hat{Y}_{0}\pm t_{\alpha/2,\nu},S_{\hat{Y}_{0}}$（公式十五）
在直角坐标系占，（公式十五）表示一条中间窄，两端宽的带子，如下所示图中的两条实曲线所示，其中最窄处对应$X_{0}=\bar{X}$，如下所示：
在R中，可以生成上述图形，由于没有用于预测的数据集，因此，只能生成一条回归直线与自变量的置信区间，如下所示：
|
|
其中划线是拟合的回归曲线，灰色的部分是自变量的置信区间。
个体Y值的预测区间
预测应时把自变量X代入到回归方程，对总体中的应变量量Y的个体值进行预测，给定X的数值$X_{0}$，对应的个体Y值也存在一个波动范围，其标准差$S_{Y_{0}}$按公式十六计算，如下所示：
$S_{Y_{0}}=S_{Y\cdot X} \sqrt{1+\frac{1}{n}+\frac{(X_{0}-\bar{X})^2}{\sum(X-\bar{X})^2}}$（公式十六）
于是给定$X=X_{0}$时，个体Y值的$1-\alpha$的预测区间为
（公式十七）
在上图中，有两条虚曲线，这两条虚曲线比实曲线范围更宽，它也是中间窄，两头宽，同样在$X_{0}=\bar{X}$处最窄。
需要注意的是，在给定$X=X_{0}$处，相应Y的均数的置信区间与其个体Y值的预测区间的含义是不同的：前者表示在固定的$X_{0}$处，如果反复抽样100次，可以算出100个相应Y的总体均数的置信区间，平均有100x(1-$\alpha$)个置信区间包含总体均数；后者表示的是一个预测值的取值范围，即预测100个个体值中平均将有100x(1-$\alpha$)个个体值在求出的范围内。
现在我们用这个案例中所求的直线回归方程，计算当$X_{0}=12$时，$\nu_{Y|X_{0}}$的95%置信区间和相应个体Y值的95%预测区间，前面我们知道，这个方程为$\hat{Y}=0.1392X + 1.667$，另外还知道$\bar{X}=9.5,l_{XX}=42$，面还计算出了$S_{Y\cdot X}=0.1970$，当$X_{0}=12$时，$\hat{Y}=1.6617+0.1392\times 12=3.3321$，根据公式十四和公式十六可得：
$S_{\hat{Y_{0}}}=0.1970\sqrt{\frac{1}{8}+\frac{(12-9.5)^2}{42}}=0.1031$
$S_{Y_{0}}=0.1970\sqrt{1+\frac{1}{8}+\frac{(12-9.5)^2}{42}}=0.2223$
查表可知$t_{0.05/2,6}=2.447$，因此按公式十五可知，当$X_{0}=12$时，尿肌酐含量总体均数的95%置信区间为：
$(3.3321-2.447\times 0.1031,3.3321+2.447\times 0.1031)=(3.080,3.584)$
按公式十七，当$X_{0}=12$时，尿肌酐含量个体值的95%预测区间为：
$(3.3321-2.447\times 0.2223,3.3321+2.447\times 0.2223)=(2.788,3.876)$
现在用R来计算一下，如下所示：
|
|
第2个案例：简单线性回归及诊断
在这个案例中，我们按照正常的直线加中归的流程计算一下，这个案例源于《医学统计学及SAS应用》（王炳顺）。
案例描述
在研究某抗心律扮演药对电刺激狗右心室致颤阈的影响，实验测得狗静脉注射不同剂量心律扮演药与右心室致颤阈的数据如下表所示：
| 药物剂量/(mg/kg) | 1 | 3 | 5 | 7 | 9 |
|---|---|---|---|---|---|
| 阈值提高/V | 8.03 | 14.97 | 19.23 | 27.83 | 36.23 |
第一步：绘制散点图
线性回归的第一步通常是绘制两个变量的散点图，在绘制的散点图中，x值表示药的剂量，y轴表示的是阈值的提高值，最终计算出来的线性方程应该是y(阈值提高值)=bx(药物剂量)+a，如下所示：
|
|
散点图如下所示：
第二步：估计回归系数
直接使用lm()函数即可，如下所示：
|
|
计算结果如下所示：
|
|
结果解释：
- 从结果中我们可以看出来，我们所求的方程y(阈值提高值)=bx(药物剂量)+a中，截矩a的值是3.9430，斜率b的值是3.4630。
- R计算的结果中还有t检验的值，与相应的p值，我们可以发现，斜率b的t检验的值小于0.001，这个值是必须要小于0.01的，也就是说必须要显著性意义，而截矩的值，不用管它。
- 线性方程的自由度是3，残差的标准差为1.448。
- $R^2$是0.9871，它表示药物的剂量可以解释98.71%的变异，后面的是调整R的平方值。
- F统计量的值为228.8，它的两个自由度为1和3，p值为0.0006272。
第三步：回归诊断
回归诊断的内容主要有三个方面：
- 误差项是否满足独立性、等方差性和正态性，选择模型是否合适。
- 是否存在异常样本，回归分析的结果是否对某些样本依赖过重，即回归模型是否具备稳定性。
- 自变量之间是否存在高度相关，即是否有多重共线性问题存在。
计算过程如下所示：
|
|
图片结果如下所示：
现在解释一下这些图形：
第一张图（左上），残差图。在构建线性回归模型时，通常采用的是最小二乘法来估计回归系数，并在此基础上作进一步的推断，前文中对其应用条件已经有强调，例如应变量与自变量的关系为线性、误差服从均数为0的正态分布，且方差相等、各观测独立等。如果实际数据不满足这些假设而直接拿来做线性回归，那么就会得到专业上无法解释的导论，至少会影响回归估计的精度与假设检验的p值，对于这些条件的检查，较为简单有效的方法就是考虑回归的残差图（residual plot）。残差图一般是将现在有模型求出来的各点残差$e_{i}=Y_{i}-\hat{Y_{i}}$作为纵坐标，相应的预测值$\hat{Y}$或者自变量取值X作为横坐标来绘制的，如果数据符合模型的基本假设，那么残差与回归预测值的散点图基本上都是呈0的上下随机分布，而不是会表现出任何特殊的结构。在这张图形中，Residual vs fitted为拟合值y对残差的图形, 可以看出,数据点都基本均匀地分布在直线y=0的两侧, 无明显趋势，满足线性假设。
现在我们看一些比较异常的残差图，如下所示：
其中a图是正常的残差图，它的残差均匀地分布在0的上面侧，说明此数据用于拟合直线回归方程是恰当的。b图，是某农药厂工人全血胆碱酯酶活性Y与工龄X进行直线回归得到的残差图，其中可以发现有一个点的残差相对于其他的点大很多，这就是一个异常点，可以考虑删除或改用其他可减小离群点影响的回归分析方法。c图是锡克反应性率Y与1-3岁儿童年龄X之间经直线回归得到的残差图，图中的残差与回归预测值表现为曲线关系，提示在目前的直线回归模型中加入自变量的二次项可能改善拟合效果（原因书中没提到，我也不明白），d图为舒张压Y与女性儿童年龄X之间直线回归的残差图，图中的残差呈喇叭口形状，说明误差的方差不剂，应考虑某种对方差进行稳定化的处理，e图表示残差之间不独立，可以看到残差与各个观测的测量时间存在较强的相关。（医学统计学.第四版.孙振球）
第二张图（右上），当预测变量值固定时，因变量成正态分布，则残差图也应是一个均值为0的正态分布。正态Q-Q图是在正态颁对应的值上，标准化残差的概率图，若满足正态假设，则图上的点应该落在45度角的直线上，若不是，则违反了正态性假设。第二幅（右上），Normal QQ-plot图中数据点分布趋于一条直线, 说明残差是服从正态分布的，用于回归的自变量的标准化残差落在(-2，2)区间以外的概率要≤0.05，若某一自变量数值的标准化残差落在(-2，2)区间以外，可在95%置信度将其判为异常实验点，不参与回归直线拟合。
第三幅图位置尺度图（Scale-Location Graph）中，水平线周围的点应随机分布，Scale-Location 图显示了标准化残差(standardized residuals)的平方根的分布情况，最高点为残差最大值点。第三副图显示基本符合方差齐性的要求。
第四幅图（Residuals vs Leverage）提供了单个观测点的信息，从图中可以鉴别离群点、高高杆值点和强影响点，主要是通过Cook距离(Cook’s distance)显示各个点对回归的影响点
还可以使用shapiro.test函数对线性回归方程的残差正态性进行验证，如下所示：
|
|
正态性检验结果表明W值为0.86452，P值为0.2449，说明残差符合正态分布。
有关残差图
在构建线性回归模型时，通常采用的是最小二乘法来估计回归系数，并在此基础上作进一步的推断，前文中对其应用条件已经有强调，例如应变量与自变量的关系为线性、误差服从均数为0的正态分布，且方差相等、各观测独立等。如果实际数据不满足这些假设而直接拿来做线性回归，那么就会得到专业上无法解释的导论，至少会影响回归估计的精度与假设检验的p值，对于这些条件的检查，较为简单有效的方法就是考虑回归的残差图（residual plot）。
残差图一般是将现在有模型求出来的各点残差$e_{i}=Y_{i}-\hat{Y_{i}}$作为纵坐标，相应的预测值$\hat{Y}$或者自变量取值X作为横坐标来绘制的，如果数据符合模型的基本假设，那么残差与回归预测值的散点图基本上都是呈0的上下随机分布，而不是会表现出任何特殊的结构。
独立性检验
判断因变量（或残差）最好的方法时依据收集数据的方式的先验知识。lmtest包提供了dwtest检验函数，car包提供了Durbin-Watson检验的函数，都能够检验误差序列的相关性，如下所示：
|
|
计算结果如下所示：
|
|
本例结果中DW值为1.9213，我不清楚这个p值的意义是什么，不过从网上的资料来看，DW值越接近于2，表示检验结果越理想。
另外一种独立性检验方法是使用car包中的crPlots()函数绘制偏残差图(partial residual plot )，也叫成分残差图(component plus residual plot)，来显示因变量与自变量之间是否呈非线性关系，如下所示：
|
|
car包中还有一个linearHypothesis函数用于线性假设检验，比图形更准，其中dose=0表示要检验自变量dose的系数是否为0，如下所示：
|
|
P值小于0.05，可以认为x的系数不为0。
方差齐性
通过以因变量为x轴，学生化残差为y轴做残差图，进行判断，如下所示：
|
|
所有的学生化残差均在$\pm 2$的范围内波动，没有明显的上升或下降趋势，可以认为符合方差齐性。还通过自变量与残差绝对值的等级相关检验来判断，如下所示：
|
|
计算结果如下所示：
|
|
自变量与残差绝对值的等级相关系数为0.2，P值大于0.05，无统计学意义，可以认为残差方差齐性。 car包提供了两个有用的函数，可判断误差方差是否恒定，ncvTest()函数生成一个计分检验，零假设为误差方差不变,备择假设为误差方差随着拟合值水平的变化而变化。若检验显著，择说明存在异方差性。spreadLevelPlot()函数创建一个添加了最佳拟合曲线的散点图，展示标准化残差绝对值与拟合值的关系，如下所示：
|
|
计分检验结果不显著，说明满足方差齐性的假设。通过水平分布图，可以看到其中的点在水平的最佳拟合曲线周围呈水平随机分布。若违反了该假设，将会呈现一个非水平的曲线。
线性模型假设的综合验证
使用gvlma中的gvlma()函数，如下所示：
|
|
计算结果如下所示：
|
|
Global Stat给模型假设提供了一个单独的综合检验（通过/不通过），本例中，可以看到数据满足线性回归模型所有的统计假设（P=0.7071）。同时还对偏度(Skewness)、峰度(Kurtosis)、连接函数(Link function)和异方差性(Heteroscedasticity)进行了评价。
影响分析
是探查对估计有异常影响的数据，如果一个样本不服从某个模型，其余数据服从这个模型，则称该样本点为强影响点。影响分析就是区分这样的样本数据。
离群点
指那些模型预测效果不佳的观测点，通常有很大的、或正或负的残差，正残差说明模型低估了响应值，负残差说明高佑了响应值。 outlierTest()函数是根据单个最大（或正或负）残差值的显著性来判断是否有离群点，若不显著，则说明数据集中没有离群点，若显著，则必须删除该离群点，然后再检验是否还有其他离群点存在。qqPlot图中落在置信区间带外的点可被认为时离群点，本例中未发现有离群点。
|
|
绘制残差的QQ图，如下所示：
|
|
高杠杆值点
高杠杆值点是与其他预测变量有关的离群点，即它们是由许多异常的预测变量组合起来的，与响应变量值没有关系。 高杠杆值的观测点可通过帽子矩阵的值（hat statistic）判断。对于一个给定的数据集，帽子均值为p/n，其中p是模型估计的参数数目（包含截距项），n是样本量。一般来说，若观测点的帽子值大于帽子均值的2或3倍，则可认定为高杠杆值点。
强影响点
强影响点，即对模型参数估计值影响有些比例失衡的点。例如，当移除 模型的一个观测点时模型会发生巨大的改变，那么需要检测一下数据中是否存在强影响点。Cook距离，或称为D统计量。Cook’s D值大于4/(n-k-1)，则表明它是强影响点，其中n为样本量大小，k是预测变量数目（有助于鉴别强影响点，但并不提供关于这些点如何影响模型的信息），如下所示：
|
|
car包中的influencePlot函数，可将离群点、杠杆点和强影响点的信息整合到一幅图形中，如下所示：
|
|
帽子统计量、DFFITS准测、Cook统计量和COVRATIO准则在R软件可分别通过hatvalues(),dffits(),cooks.distance()和covration()函数计算。influence.measures()可对一次获得这四个统计量的结果，如下所示：
|
|
纵坐标超过2或小于-2的值可被认为是离群点，水平轴超过0.2或0.3的值有高杠杆值（通常为预测值的组合）。圆圈大小与影响成比例，圆圈很大的点可能是对模型估计造成的不成比例影响的强影响点，influence.measures()的inf用星号标注异常值。
多元线性回归
直线回归是研究一个应变量与单个自变量之间呈直线关系的一种统计方法。但是在医学研究中，许多疾病都有多种原因，而且预后也是多种因素决定的，如糖尿病患者的血糖变化可能受胰岛素、糖化血红蛋白、血清总胆固醇、甘油三酯等多种生化指标的影响，这个时候就需要研究多个因素对一个结果变量数量依存的线性关系以及多个因素之间的线性关系，这种关系就称为多元线性回归。例如如人体的体表现积随身高和体重而改变，则体表面积为应变量，身高和体重为自变量，设体表现为Y，身高和体重为X1和X2，则所建立的含有两个自变量二元回归方程为：
$\hat{y}=a+b_{1}x_{1}=b_{2}x_{2}$
其中$b_{1}\hat{y}$是y均数的估计值或预测值（predicted value，y hat），a为截距，表示各自变量取0时，y的估计值，$b_{1}b_{2}$和$b_{2}b_{1}$称为偏回归系数（partial regression coefficient），其中$b_{1}$表示在体重的取值固定不变时，身高每增加一个单位，y估计值平均变化$b_{1}$个单位；$b{_2}$表示在身高的取值不变时，体重每增加一个单位，y估计值平均变化$b_{2}$个单位。一般来讲，对于拟合一个应变量y和m个自变量$x_{1},x_{2},\cdots,x_{m}$的多元回归方程，其形式为：
$\hat{y}=a+b_{1}x_{1}+b_{2}x_{2}+\cdots=b_{i}x_{i}+\cdots+b_{m}x_{m}$
其中，a为截距，意义同上，$b_{i}$为当其他自变量取值固定不变时，$x_{i}$每增加一个单位，y的估计值平均变化$b_{i}$个单位。
多元线性回归资料应满足以下条件：
- 自变量和应变量之间的关系是线性关系；
- 各观测单位相互独立；
- 残差服从正态分布；
- 残差满足方差齐性。
多元线性回归的原理
对于多元回归方程的建立，也就是求出各常数项和各偏回归系数，原理还是使用最小二乘法来求各待估计系数，即根据观察到的n例数据，使残差平方和最小，即如下所示：
$Q=\sum_{k=1}^n(Y_{k}-\hat{Y_{k})^2}=\sum_{k=1}^n[Y_{k}-(b_{0}+b_{1}X_{1k}+b_{2}X_{2k}+\cdots=b_{m}X_{mk})]^2$（公式十八）
由此可以得到下面的方程组，求解得到$b_{1},b_{2},\cdots,b_{m}$,并求出$b_{0}$项，如下所示：
$b_{0}$的计算如下所示：
这两个公式分别表示自变量的离均差平方和（i=j），两个自变量的离均差积和（$i\neq j$）及自变量$X_{j}$与应变量Y的离均差和和。
多元线性回归的计算过程
现在看一个案例：
27名糖尿病患者的血清总胆固醇、甘油三酯、空腹胰岛素、糖化血红蛋白、空腹血糖的测量值位于下表中，试建立血糖与其他几项指标的多元线性回归方程（《医学统计学》第四版.孙振球，第15章）。
第一步，由上表的数据，通过公式（二十一）和公式（二十二）计算出包括应变量在内的各变量离差矩阵，如下所示：
根据公式（十九）列出方程组求解后，可以得到：
b1=0.1424, b2=0.3515, b3=-0.2706, b4=0.6382
算得均数分别为：
$\bar{X_{1}}=5.8126,\quad \bar{X_{2}}=2.8407,\quad\bar{X_{3}}=6.1467,\quad\bar{X_{4}}=9.1185,\quad \bar{Y}=11.9259$
根据公式（二十）可以计算出截矩，如下所示：
$b_{0}=11.9259-(0.1424 \times 5.8126+0.3515 \times 2.8407-0.2706 \times6.1467 +0.6382 \times9.1185)=5.9433$
有了这几个系数与截矩，我们就可以得到线性回归方程，如下所示：
$\hat{Y}=5.9433+0.1424X_{1}+0.3515X_{2}+ -0.2706X_{3}+0.6382X_{4}$
在R中计算多元线性回归方程与简单线性回归基本上类似，也是使用lm()函数，如下所示：
|
|
计算结果如下所示：
|
|
从结果中可以看出来，截矩，X1，X2，X3，X4的偏回归系数分别为5.9433，0.1424，0.3515，-0.2706，0.6382，对回归的结果进行方差分析可以看出来，回归方程的系数的显著性不高, 有的没有通过检验，这说明如果选择全部变量构造方程, 效果并不好，这就涉及到变量选择的问题，以建立“最优”的回归方程。
如果先不考虑变量的筛选问题，那么我们计算得到的多元线性回归方程就是下面的这个样子：
$\hat{Y}=5.9433+0.1424X_{1}+0.3515X_{2}-0.2706X_{3}+0.6382X_{4}$
多元线性回归中变量的筛选
在R中，进行“最优”回归方程的方法是“逐步回归法”的计算函数step( ),它是以Akaike信息统计量为准则(简称AIC准则), 通过选择最小的AIC信息统计量, 来达到删除或增加变量的目的。AIC的数值越小，就说明对变量的筛选效果越好，它的用法如下所示：
|
|
其中，object是线性模型或广义线性模型的结果，scope是确定逐步回归的搜索区域，direction确定逐步搜索的方向，both是一切子集回归法，backward是向后法，forward是向前法，默认值为both，现在使用setp()对上述的多元线性回归结果做逐步回归，如下所示：
|
|
计算结果如下所示：
|
|
结论: 用全部变量作回归方程时, AIC统计量的值为42.16,如果去掉变量X1，AIC统计量的值为40.34，如下去掉变量X2，AIC统计量的值为43.568，依次类推。由于去掉X1命名AIC统计量达到最小，因此R会自动掉变量X1，进入下一轮的计算，在下一轮中，无论去掉哪个变量，AIC的统计量的值均会升高，因此R会自动终止计算，得到“最优”回归方程，现在使用summary()提取相关回归信息，如下所示：
|
|
结论：从上面的结果我们可以看出来，回归系数的显著性水平有很大提高，每个自变量的检验均是显著的，𢅟得到了“最优”的回归方程，如下所示：
$\hat{Y}=6.4996+0.4023X_{2}-0.2870X_{3}+0.6632X_{4}$
回归诊断
由样本数据得到回归方程后，为了确定回归方程及引入的自变量是否有统计学意义，必须进一步做假设检验。方差分析法可以将回归方程中所有自变量$X_{1},X_{2},\cdots,X_{m}$作为一个整体来检验它们与应变量Y与之间是否具有线性关系，并对回归方程的预测或解释能力作出综合评价。在此基因上进一步对各变量的重要性作出评价。
方差分析
方差分析法原理
在多元线性回归中，为什么要做方差分析呢，这是因为在线性回归中，把因变量的总变异，也就是因变量的所有观测值到其均值$\bar{y}$的距离平方和（通常用SST）表示，分解成对应于各个自变量的若干平方和(SS,sum of squares)，每一个这样的平方和代表了一个自变量对总变异平方的贡献，剩下不能被自变量解释的，就是误差平方和或残差平方和(SSE,error sum of squares)，自变量及残差平方和所对应的平方和加起来等于因变量的总变异。
用简单的话来讲，就是方差分析表把由自变量所导致的因㸄的变化和误差所导致的因变量的变化找出来，并进行比较（通过F检验），如果这两个变化很不一样，那么则认为自变量在回归中的确显著影响了因变量，表示回归有道理，否则，没有理由认为因变量和误差有所不同。
具体来讲，如果我们有两个自变量，分别为A和B，则把SSt分解为对应于这两个自变量的两个平方和SSA，SSB以及残差平方和SSE，即SST=SSA+SSB+SSE。那么如果SSA与SSE相比较很显著，即变量A的贡献显著大于误差的贡献，则称变量A显著，否则，如果SSA和SSE的贡献差不多，那么变量A就不显著。这种显著需要正式的检验，如果前面关于误差的四个假定成立（回归模型的误差项独立，且服从正态分布），那么这些平方和（比如这里的SSA、SSB、SSE）都有$\chi^2$分布（因为它们是正态分布变量的平方和），而除以各自的自由度之后，得到它们的平均（叫均方，mean square，记为MSA，MSB，MSE），根据F分布的定义，MSA/MSE和MSB/MSE都有F分布，如果MSA/MSE很大，那么A就是和残差相比显著，则相应的p值就小。
方差分析法过程
方程的假设检验如下所示：
$H_{0}:\beta_{1}=\beta_{2}=\cdots=\beta_{m}=0$
$H_{1}$各$\beta_{j}(j=1,2,\cdots,m)$不全为0
将应变量Y的总变异分解成两部分，即如下所示：
$\sum(Y-\bar{Y})^2=\sum(\hat{Y}-\bar{Y})+\sum(Y-\hat{Y})^2$
其中$\sum(Y-\bar{Y})^2$ 是回归平方和，$\sum(Y-\hat{Y})^2是残差平方和，上式可记作：
$SS_{总}=SS_{回}+SS_{残}$
其中，回归平方和可用下式计算：
$SS_{回}=b_{1}l_{1Y}+b_{1}l_{2Y}+\cdots+b_{m}l_{mY}=\sum b_{j}l_{jY}$
残差平方和为：$SS_{残}=SS_{总}-SS_{回}$
用统计量F检验假设$H_{0}$是否成立，如下所示：
$F=\frac{SS_{回}/m}{SS_{残}/(n-m-1)}=\frac{MS_{回}}{MS_{残}}$
方差分析见下表：
| 变异来源 | 自由度 | SS | MS | F |
|---|---|---|---|---|
| 总变异 | n-1 | SS总 | ||
| 回归 | m | SS回 | SS回/m | MS回/MS残 |
| 残差 | n-m-1 | SS残 | SS残/(n-m-1) |
如果$F\geq F_{\alpha,(m,n-m-1)}$，则在α水准上拒绝H0，接受H1，认为应变量Y与m个自变量$X_{1},X_{2},\cdots,X_{m}$之间存在线性回归关系。现在我们从刚才的多元线性回归方程中计算各部分的变异：
SS总=222.5519
SS回=0.1424×67.6962+0.3515×89.8025+0.2706×142.4337+0.6382×84.5570=133.7107
SS残=222.5519-133.7107
方差分析的结果如下表所示：
| 变异来源 | 自由度 | SS | MS | F | P |
|---|---|---|---|---|---|
| 总变异 | 26 | 222.5519 | |||
| 回归 | 3 | 133.0978 | 44.36593 | 11.4 | <0.01 |
| 残差 | 23 | 89.4541 | 3.89 |
查F界值表可知，$F_{0.01,(3,23)}=4.76$，F>4.76，P<0.01，在α=0.05检验水准上，拒绝H0，接受H1，认为所拟合的回归方程具有统计学意义。
注，在R中查找相应F值的函数为qf()，用法如下所示：
|
|
使用我要查找在0.01水平上，两个自由度为3和23的F值，如下所示：
|
|
在R中，可以通过anova()来查看多元线性回归的方差分析结果，如下所示：
|
|
决定系数$R^2$
根据方差分析的结果，还可以得到多元线性回归的决定系数$R^2$，它的计算公式如下所示：
$R^2=\frac{SS_{回}}{SS_{总}}=1-\frac{SS_{残}}{SS_{总}}$
从上面的结果可以看出来$R^2=0.5981$，这个数字从前面的方差分析结果中也能看出来，即Multiple R-squared: 0.5981。这个数字表明，血糖含量变量的58%可由甘油三脂(X2)、胰岛素(X3)、糖化血红蛋白(X4)的变化来解释。
复相关系数
$R=\sqrt{R^2}$称为复相关系数（multiple correlation coefficient），可用来度量应变量Y与多个自变量的线性相关程度，亦即观察值Y与估计值$\hat{Y}$之间的相关程度，在这个案例中，复相关系数$R=\sqrt{0.5981}=0.7734$，如果只有一个自变量时，$R=|r|$，r为简单相关系数。
t检验
方差分析和决定系数是将所有自变量$X_{1},X_{2},\cdots,X_{m}$作为一个整体来检验和说明它们与Y的相关程度及解释能力的，并未指明方程中的每一个自变量对Y的影响如何，但是在医学研究中，我们往往更加关心是对各自变量的解释。对每一自变量的作用进行检验和称量它们对Y的作用大小，可以采用偏回归平方和法，t检验法，标准化回归系数法。
而在R中的，通常采用的是t检验法，在前面的多元线性回归中，我们就可以看到这个t检验法的结果，如下所示：
|
|
根据t界值表，我们可以发现，剔除了X1变量后，剩下的变量的p值都小于0.05（在R中，一个星号表示大于0.01，小于0.05）,说明这些变量有统计学意义，对于同一资料而言，不同自变量的t值间可以相互比较，t的绝对值越大，说明该自变量对Y的回归所起的作用越大。
残差分析
残差分析是检查资料是否符合模型条件的一个有用的工具。如果资料与模型的假设偏离较大，最小二乘估计、假设检验及区间估计则有可能出现问题。绘制多元线性回归的残差图，如下所示：
|
|
共线性诊断
共线性问题是指拟合多元线性回归时，自变量之间存在线性关系或氛这线性关系，自变量之间的线性关系将会隐藏变量的显著性，增加参数估计的误差，还会产生一个很不稳定的模型，因此共线性诊断就是要找出哪些变量之间存在共线性关系，共线性诊断的方法主要有特征值法，条件指数法，方差膨胀因子法。具体的原理可以参考相应的统计学书籍，例如《R语言与统计分析》（汤银才），这里只说一种方法，即方差膨胀因子法，使用的函数为car包中的vif()，VIF的全称是Variance Inflation Factor，即方差膨胀因子，一般原则下，如果VIF值大于4就表明存在多重共线性问题，它的使用如下所示：
|
|
从结果我们可以看出，这几个变量的VIF值都没有大于4，因此不存在多重共线性问题。
参考资料
- 孙振球. 医学统计学.第4版[M]. 人民卫生出版社, 2014.
- DawnGriffiths. 深入浅出统计学[M]. 电子工业出版社, 2012.
- 王炳顺. 医学统计学及SAS应用[J]. 2007.
- 王斌会. 多元统计分析及R语言建模[M]. 暨南大学出版社, 2010.
- 汤银才. R语言与统计分析[M]. 高等教育出版社, 2008.
- 卡巴科弗. R语言实战[M]. 人民邮电出版社, 2016.