多元線性回歸分析示例

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
           

運作結果如圖所示：

結果分析：