Python和R中多元线性回归的实证分析

回归分析方法说白了就是处理多个变量相互依赖关系的一种数理统计方法(之前并没学过数理统计，恶补了一下，挺爽的~)。这篇随笔中主要运用了线性代数和数理统计知识，欢迎各方大佬指正,错误之处，不胜感激。

一.建立模型

这里我们假定研究变量Y与x₁,x₂,x₃……x_m，m个变量之间的相互依赖关系。采取现实生活中观测的n组变量Y与变量x数据，建立如下方程组：

y_i=β0+β₁x_i1+β₂x_i2+……β_mx_im+ε_i(i=1,2,3……n)

即：Y=Cβ+ε

为弥补建立的方程组与实际数据的误差，引入ε为随机的误差变量，服从均值为0，方差为1的标准正态分布。Y为n个y观测值得列向量，C为n行m+1列的x观测值矩阵，ε为n个误差变量的列向量。

二.模型完善

1.β的最小二乘估计

将上述模型中的β的最小二乘估计量(最小二乘法：数理统计中利用微分知识，求得误差平方和最小的方法)设为b=(b₀,b_1,……,b_m)^'.

此时Q(b)=minQ(β)(对于一切β),其中Q(b)=∑εi²=∑[yi-(β0+β₁x_i1+β₂x_i2+……β_mx_im+ε_i)]²=(Y-Cβ)^'(Y-Cβ).

由数学推导(硬伤，这里直接拿来用了)可得β的最小二乘估计量β^{^}=b=(C^'C)^-1C^'Y(方程有解的前提下，C的秩为m+1)。设Y^{^}=Cβ^{^}=C(C^'C)^-1C^'Y=HY为Y的预测向量，其中H_n×n=C(C^'C)^-1C^'被称为奇怪的帽子矩阵(无法理解数学家们的思维，长得像帽子- -!)。此时残差向量为

ε^{^}=Y-Y^{^}

且残差平方和为

Q(β^{^})=ε^{^}'ε^{^}

2.σ2的估计

　利用最大似然比原理(拿来用就好),可得β的最大似然估计量仍为b，　σ²的最大似然估计为

σ^{^2}=1/nΣ[yi-(b0+b₁x_i1+b₂x_i2+……b_mx_im)]²

　  =1/n(Y-Cb)^'(Y-Cb)=1/nQ(B)

但σ^{^2}不是σ²的无偏估计量(取得的样本期望与实际真值相同)，取

s²=1/(n-m-1)Q(b)......(样本方差)

作为σ²的估计量(无偏估计量)。

3.回归方程的显著性检验

建立模型时我们假设了Y与m个自变量之间存在线性关系。模型只是一个假设，为此在求解出各无偏估计量之后要对Y与自变量之间是否存在线性关系进行检验。

实际上，该问题可化为对求解出的β矩阵中的各系数进行检验，若Y与m个自变量之间均无线性关系，则β为全零矩阵。为选择合适检验统计量，引入平方和分解公式(原理很简单)。

(1)平方和分解公式

给定观测矩阵：

恒有平方和分解公式公式：

Σ(y_i-y)2=Σ(y_i-y_i^)2+Σ(y_i2-y)2

其中

y=1/nΣy_i

y_i^=Cβ^{^}

平方和分解公式公式左边体现了Y的观测值总波动大小，称为总偏差平方和，记作l_yy(或TSS)。平方和分解公式公式右侧第二项体现了n个无偏估计量的波动大小，它是由于Y与各自变量之间存在线性关系，并通过自变量的变换而引起称为回归平方和，记作U(或MSS)；平方和分解公式公式右侧第一项称作残差平方和(误差平方和)，记作Q(或ESS)。

　基于上述说明，平方和分解公式公式可简记为：

  l_yy=Q+U或TSS=ESS+MSS

(2)相关性检验

由平方和分解公式及MSS和ESS的意义，若MSS比ESS大得多，则Y的总偏差TSS主要由自变量的变化引起的。比值MSS/ESS可以作为检验统计量。

构造如下检验统计量为：

F = (MSS/m)/(ESS/(n-m-1))=模型均方/均方误差。(我猜这里各自比上自己的*的作用是归一化处理)

三.模型求解

多元线性回归在Matlab中主要实现函数如下(自己查看Matlab中的用户手册吧，最强大的老师)：

(1) b=regress(Y,X)确定回归系数的点估计值

(2) [b,bint,r,rint,stats]=regress(Y,X,alpha)

(3) rcoplot(r,rint)画出残差及其置信区间

实例：为了分析各类人群血压的差异及其原因，医学部门经过大量调查，认为血压主要与年龄、体重指数、吸烟习惯等因素有关。根据已给数据，建立血压与各因素之间的回归模型。

分析：建立Y矩阵,自变量X矩阵，调用Matlab库函数。贴代码

n=30;m=3;
Y=[144 215 138 145 162 142 170 124 158 154 162 150 140 110 128 130 135 114 116 124 136 142 120 120 160 158 144 130 125 175];%27个人各自的血压值
x1=[39 47 45 47 65 46 67 42 67 56 64 56 59 34 42 48 45 18 20 19 36 50 39 21 44 53 63 29 25 69];%27个人各自的年龄
x2=[24.2 31.1 22.6 24.0 25.9 25.1 29.5 19.7 27.2 19.3 28.0 25.8 27.3 20.1 21.7 22.2 27.4 18.8 22.6 21.5 25.0 26.2 23.5 20.3 27.1 28.6 28.3 22.0 25.3 27.4];%27个人各自的体重指数
x3=[0 1 0 1 1 0 1 0 1 0 1 0 0 0 0 1 0 0 0 0 0 1 0 0 1 1 0 1 0 1];%27个人各自的吸烟习惯
X=[ones(n,1),x1',x2',x3'];
[b,bint,r,rint,s]=regress(Y',X);%b为回归系数，bint为回归系数的置信区间，s有三个值，分别是相关%%系数，检验统计量F，以及与F对应的概率
s2=sum(r.^2)/(n-m-1);%样本方差(方差无偏估计量)
b,bint,s,s2
rcoplot(r,rint)

 运行结果：

b =

   45.3636     0.3604     3.0906    11.8246

bint =

    3.5537   87.1736    -0.0758    0.7965     1.0530    5.1281    -0.1482   23.7973

s =

    0.6855   18.8906    0.0000  169.7917

s2 =

  169.7917

笔者实在贴不上图片，残差图自行脑补。由残差图上颜色可知第二、十个点残差过大，将其剔除后，运行结果如下：

b =

   58.5101     0.4303     2.3449    10.3065

bint =

   29.9064   87.1138     0.1273    0.7332     0.8509    3.8389     3.3878   17.2253

s =

    0.8462   44.0087    0.0000   53.6604

s2 =

   53.6604

方差大幅度降低，由此可建立血压与各因素之间的回归模型为：

y=58.5101+0.4303x₁+2.3449x₂+10.3065x₃.