GLM模型应用于脑功能影像分析时，在某个因素影响下，由beta图，经过t检验得到脑区显著激活的区域。应用于其他地方也可加深我们对于模型的理解。
clc,clear;
X=[ 136.5 215136.5 250136.5 180138.5 250138.5 180138.5 215138.5 215138.5 215140.5 180140.5 215140.5 250];
y=[ 6.27.54.85.14.64.64.94.12.83.14.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<perror('观察值的数目过少');
end
if nargin < 2error('多元线性回归要求有两个输入参数');
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~=pdisp('预测自变量的个数不正确');returnendend
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}='显著';
elsetable{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:ptable(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
运行结果如图所示：
结果分析：
