多元线性回归分析示例

GLM模型应用于脑功能影像分析时，在某个因素影响下，由beta图，经过t检验得到脑区显著**的区域。应用于其他地方也可加深我们对于模型的理解。

clc,clear;
X=[     136.5          215
        136.5          250
        136.5          180
        138.5          250
        138.5          180
        138.5          215
        138.5          215
        138.5          215
        140.5          180
        140.5          215
        140.5          250];
y=[       6.2
          7.5
          4.8
          5.1
          4.6
          4.6
          4.9
          4.1
          2.8
          3.1
          4.3 ];
Xnew=[137.5,240];
pp=0.95;
[ab,stats,yy,ylr]=regres2(X,y,Xnew,pp)
table=stats{1}

调用的回归函数如下 ；

function [beta,stats,ynew,ylr]=regres2(X,y,Xnew,pp)
beta=[];stats=[];ynew=[];ylr=[];
[n,p]=size(X);m=p+1;
if n<p
    error('观察值的数目过少');
end
if  nargin < 2
    error('多元线性回归要求有两个输入参数');
end 
[n1,collhs] = size(y);
if n ~= n1, 
    error('输入参数y的行数，必须等于输入参数X的行数.'); 
end 
if collhs ~= 1, 
    error('输入参数y应该是一个列向量'); 
end
if nargin==3   
    if isnumeric(Xnew)
        [n1,p1]=size(Xnew);
        if p1~=p
            disp('预测自变量的个数不正确');
            return
        end
    end
end
if (nargin<4)|(~isnumeric(pp))|(pp<=0)|(pp>=1)
    pp=0;
end
A=[ones(size(y)),X];
[beta,btm1,rtm,rtm1,stat] =regress(y,A);
alpha=[0.05,0.01];
yhat=A*beta;
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=finv(1-alpha,m-1,n-m);
table=cell(p+4,7);
table(1,:)={'方差来源','偏差平方和','*度','方差','F比','Fα','显著性'};
table(2+p,1:6)={'回归',SSR,m-1,SSR/(m-1),Fb,min(Falpha)};
table(3+p,1:6)={'剩余',SSE,n-m,SSE/(n-m),[],max(Falpha)};
table(4+p,1:3)={'总和',SST,n-1};
if Fb>max(Falpha)
    table{2+p,7}='高度显著';
elseif (Fb<=max(Falpha))&(Fb>min(Falpha))
    table{2+p,7}='显著';
else
    table{2+p,7}='不显著';
end
R2=SSR/SST;R=sqrt(R2);
Sy=sqrt(SSE/(n-m));
mnX=mean(X);
MNX=repmat(mnX,n,1);
Ljj=diag((X-MNX)'*(X-MNX));
Pj=abs(beta(2:end).*sqrt(Ljj/SST));
C=diag(inv(A'*A));bj2=beta.*beta;
SSj=bj2(2:end)./C(2:end);
Fj=SSj/SSE*(n-m);
Falpha=finv(1-[0.05,0.01],1,n-m);
ind2=find(Fj>=Falpha(2));
ind1=find((Fj>=Falpha(1))&(Fj<Falpha(2)));
ind0=find(Fj<Falpha(1));
xxx=zeros(size(Fj));
xxx(ind2)=2;
xxx(ind1)=1;
[tmp,zbx]=min(Fj);
xzh={'不显著','显著','高度显著'};
for kk=1:p
    table(kk+1,:)={['x',num2str(kk)],SSj(kk),1,SSj(kk),Fj(kk),[],xzh{1+xxx(kk)}};
end
table{2,6}=Falpha(1);table{3,6}=Falpha(2);
stats={table,R,Sy,Pj};
if (nargin>2)&(isnumeric(Xnew))
    [n1,p1]=size(Xnew);
    Xnew=[ones(n1,1),Xnew];
    ynew=Xnew*beta;
    Shat2=SSE/(n-m)*(1+Xnew*inv(A'*A)*Xnew');
    Syhat=sqrt(diag(Shat2));
    ta=tinv(0.5+pp/2, n-p-1);
    yl=ynew-ta*Syhat;
    yr=ynew+ta.*Syhat;
    ylr=[yl(:),yr(:)];
end

运行结果如图所示：

结果分析：