多元线性回归数据筛选问题

前几天，偶然在网上看到，建立多元线性回归时须对其中的异常数据进行筛选剔除，那这里的异常点指的是什么呢？这里的异常点指的是人为采集数据误差或者某些异常的个例等等一些不太准确的数据。

例子说明

例如，对于下表中的一些数据我们来进行数据筛选。
在医学上，糖尿病人的血糖量y(mmol/L)$y(mmol/L)$与总胆固醇 1 X (mmol / L)，甘油三脂X2(mmol/L)$X_2(mmol/L)$胰岛素X3(μU/ml)$X_3(\mu U/ml)$糖化血红蛋白X4($X_4($等有关。附表中是某医院的原始数据。我们将对其中的异常数据进行筛选。

血糖	总胆固醇	甘油三脂	甘油三脂	胰岛素
11.2	5.68	1.9	4.53	8.2
8.8	3.79	1.64	7.32	6.9
12.3	6.02	3.56	6.95	10.8
11.6	4.85	1.07	5.88	8.3
13.4	4.6	2.32	4.05	7.5
18.3	6.05	0.64	1.42	13.6
11.1	4.9	8.5	12.6	8.5
12.1	7.08	3	6.75	11.5
9.6	3.85	2.11	16.28	7.9
8.4	4.65	0.63	6.59	7.1
9.3	4.59	1.97	3.61	8.7
10.6	4.29	1.97	6.61	7.8
8.4	7.97	1.93	7.57	9.9
9.6	6.19	1.18	1.42	6.9
10.9	6.13	2.06	10.35	10.5
10.1	5.71	1.78	8.53	8
14.8	6.4	2.4	4.53	10.3
9.1	6.06	3.67	12.79	7.1
10.8	5.09	1.03	2.53	8.9
10.2	6.13	1.71	5.28	9.9
13.6	5.78	3.36	2.96	8
14.9	5.43	1.13	4.31	11.3
16	6.5	6.21	3.47	12.3
13.2	7.98	7.92	3.37	9.8
20	11.54	10.89	1.2	10.5
13.3	5.84	0.92	8.61	6.4

具体做法如下：
我们先通过Matlab$Matlab$对原始数据进行检验，对残差进行分析，得到了残差分析图，剔除其中的异常点。

从图上可以看出，第13$13$个点和第27$27$个点是异常点，这样在做数据处理时就可以将其剔除。
下面是我遇到的一些问题：

• 至今没有找到这个算法的数学原理，没看到书上有残差向量的相关介绍。
• Matlab$Matlab$程序是改进网上的开源程序而来的，从程序中看不出他为什么要这么做。
• 这种筛选数据的方法是否只适应于线性回归，非线性回归是否也可以借鉴。

Matlab代码

clear all;
clc;

%输入数据
z=xlsread('data.xls');
z1=z;
y=z(:,1);
X=[ones(size(y)),z(:,2:5)];
alpha=0.05;
b=inv((X'*X))*X'*y;  %回归系数
[n,ncolX] = size(X);

%剔除异常数据
wasnan=(isnan(y)|any(isnan(X),2));
havenans=any(wasnan);
[Q,R,perm]=qr(X,0);
p=ncolX;
RI=R\eye(p);
nu=max(0,n-p);
yhat=X*b;
r=y-yhat;  %残差
normr=norm(r);
rmse=normr/sqrt(nu);    % Root mean square error.
tval=tinv((1-alpha/2),nu);
s2=rmse^2;
hatdiag=sum(abs(Q).^2,2);
ok=((1-hatdiag)>sqrt(eps(class(hatdiag))));
denom=(nu-1).*(1-hatdiag);
sigmai=zeros(length(denom),1);
sigmai(ok)=sqrt(max(0,(nu*s2/(nu-1))-(r(ok).^2./denom(ok))));
ser=sqrt(1-hatdiag).*sigmai;
rint=[(r-tval*ser) (r+tval*ser)];
rcoplot(r,rint);
kk=[];
for i=1:n
    if (rint(i,1)>0&&rint(i,2)>0)||(rint(i,1)<0&&rint(i,2)<0)
       kk=[kk,i]; 
    end
end
X(kk,:)=[];
y(kk)=[];

%剔除异常点后，求解回归系数
beta=inv((X'*X))*X'*y;  %回归系数

%回归分析
X1=X;
X1(:,1)=[];
n=size(y,1);    %观察单位数
m=size(X,2);    
p1=m-1;          %自变量个数

alpha=0.05;  
yhat=X*beta;

%方差分析表（F检验）
SSR=(yhat-mean(y))'*(yhat-mean(y));  %回归平方和

SSE=(yhat-y)'*(yhat-y);        %残差平方和

SST=(y-mean(y))'*(y-mean(y));    %总平方和

Fb=(SSR/(m-1))/(SSE/(n-m));       %显著性检验的统计量

Falpha=2*(1-fcdf(abs(Fb),m-1,n-m));

table1=cell(4,6);    %创建元胞
table1(1,:)={'模型','偏差平方和','*度','均方','F值','F.Sig'};
table1(2,1:6)={'回归',SSR,m-1,SSR/(m-1),Fb,Falpha};
table1(3,1:6)={'残差',SSE,n-m,SSE/(n-m),[],[]};
table1(4,1:3)={'总和',SST,n-1};

C=diag(inv((X'*X)));
bj2=beta.*beta;     %回归系数平方
SSj=bj2(2:end)./C(2:end);     %偏回归系数平方和

%决定系数检验
R2=SSR/SST;   %决定系数
R=sqrt(R2);   %复相关系数
RC=1-(1-R2)*(n-1)/(n-1-p1);
Sy=sqrt(SSE/(n-m));   %剩余标准差
table2=cell(2,5);    %创建元胞
table2(1,:)={'模型','R','R平方','校正决定系数','剩余标准差'};
table2(2,1)={1};
table2(2,2)={R};
table2(2,3)={R2};
table2(2,4)={RC};
table2(2,5)={Sy};

%t检验
s=zeros(m,1);       %回归系数的标准误差
for i=1:m
    s(i,1)=sqrt(C(i))*sqrt(SSE/(n-m));
end

mnX=mean(X1);
MNX=repmat(mnX,n,1);   %复制mnX到矩阵MNX中
Ljj=diag((X1-MNX)'*(X1-MNX));     %Ljj的对角线元素为(X1-MNX)'*(X1-MNX)
Pj=beta(2:end).*sqrt(Ljj/SST);   %标准偏回归系数；

t=zeros(m,1);
for i=1:m
    t(i,1)=beta(i,1)/s(i,1);
end

p2=zeros(m,1);
for i=1:m
    p2(i,1)=2*(1-tcdf(abs(t(i,1)),n-m));
end

table3=cell(m+1,6);    %创建元胞
table3(1,:)={'模型','偏回归系数','回归系数的标准误差','标准偏回归系数','t值','P值'};
table3(2,:)={'常量',beta(1,1),s(1,1),[],t(1,1),p2(1,1)};

for kk=1:m-1
    table3(kk+2,:)={['x',num2str(kk)],beta(kk+1,1),s(kk+1,1),Pj(kk,1),t(kk+1,1),p2(kk+1,1)};
end

disp('系数分析表');
disp(table3);
disp('模型汇总');
disp(table2);
disp('方差分析表');
disp(table1);