SurfStatF

F statistics for comparing two uni- or multi-variate fixed effects models.  

 Usage: slm = SurfStatF( slm1, slm2 );

 slm1 and slm2 must have the following fields:
 slm.X   = n x p design matrix.
 slm.df  = degrees of freedom.
 slm.SSE = k*(k+1)/2 x v matrix of sum of squares of errors.

 All fields of the bigger model are copied to slm, and 
 slm.k  = k=#variates.
 slm.df = [df1-df2, df2], where df1 and df2 are the min and max of slm1.df
          and slm2.df, and SSE1 and SSE2 are the corresponding slm_.SSE's.
 slm.t  = l x v matrix of non-zero eigen values, in descending order, of 
          F = (SSE1-SSE2)/(df1-df2)/(SSE2/df2), where l=min(k,df1-df2);  
          slm.t(1,:) = Roy's maximum root = maximum F over all linear 
          combinations of the k variates. k>3 is not programmed yet.

0001 function slm = SurfStatF( slm1, slm2 );
0002 
0003 %F statistics for comparing two uni- or multi-variate fixed effects models.
0004 %
0005 % Usage: slm = SurfStatF( slm1, slm2 );
0006 %
0007 % slm1 and slm2 must have the following fields:
0008 % slm.X   = n x p design matrix.
0009 % slm.df  = degrees of freedom.
0010 % slm.SSE = k*(k+1)/2 x v matrix of sum of squares of errors.
0011 %
0012 % All fields of the bigger model are copied to slm, and
0013 % slm.k  = k=#variates.
0014 % slm.df = [df1-df2, df2], where df1 and df2 are the min and max of slm1.df
0015 %          and slm2.df, and SSE1 and SSE2 are the corresponding slm_.SSE's.
0016 % slm.t  = l x v matrix of non-zero eigen values, in descending order, of
0017 %          F = (SSE1-SSE2)/(df1-df2)/(SSE2/df2), where l=min(k,df1-df2);
0018 %          slm.t(1,:) = Roy's maximum root = maximum F over all linear
0019 %          combinations of the k variates. k>3 is not programmed yet.
0020 
0021 if isfield(slm1,'r') | isfield(slm2,'r')
0022     warning('Mixed effects models not programmed yet.');
0023     return
0024 end
0025 
0026 if slm1.df>slm2.df
0027     X1=slm1.X;
0028     X2=slm2.X;
0029 else
0030     X1=slm2.X;
0031     X2=slm1.X;
0032 end
0033 
0034 r=X1-X2*pinv(X2)*X1;
0035 d=sum(r(:).^2)/sum(X1(:).^2);
0036 if d>eps
0037     error('Models are not nested')
0038     return
0039 end
0040 
0041 if slm1.df>slm2.df
0042     df1=slm1.df(length(slm1.df));
0043     df2=slm2.df(length(slm2.df));
0044     SSE1=slm1.SSE;
0045     SSE2=slm2.SSE;
0046     slm=slm2;
0047 else
0048     df1=slm2.df(length(slm2.df));
0049     df2=slm1.df(length(slm1.df));
0050     SSE1=slm2.SSE;
0051     SSE2=slm1.SSE;
0052     slm=slm1;
0053 end
0054 
0055 slm.df=[df1-df2, df2];
0056 h=SSE1-SSE2;
0057 if ndims(slm.coef)==2
0058     slm.k=1;
0059     slm.t=h./(SSE2+(SSE2<=0)).*(SSE2>0)*(df2/(df1-df2));
0060 else
0061     [k2,v]=size(SSE2);
0062     k=round((sqrt(1+8*k2)-1)/2);
0063     slm.k=k;
0064     if k>3
0065         warning('Roy''s max root for k>3 not programmed yet.');
0066         return
0067     end
0068     l=min(k,df1-df2);
0069     slm.t=zeros(l,v);
0070     if k==2
0071         det=SSE2(1,:).*SSE2(3,:)-SSE2(2,:).^2;
0072         a11=SSE2(3,:).*h(1,:)-SSE2(2,:).*h(2,:);
0073         a21=SSE2(1,:).*h(2,:)-SSE2(2,:).*h(1,:);
0074         a12=SSE2(3,:).*h(2,:)-SSE2(2,:).*h(3,:);
0075         a22=SSE2(1,:).*h(3,:)-SSE2(2,:).*h(2,:);
0076         a0=a11.*a22-a12.*a21;
0077         a1=(a11+a22)/2;
0078         s1=real(sqrt(a1.^2-a0));
0079         d=(df2/(df1-df2))./(det+(det<=0)).*(det>0);
0080         slm.t(1,:)=(a1+s1).*d;
0081         if l==2
0082             slm.t(2,:)=(a1-s1).*d;
0083         end
0084     end
0085     if k==3
0086         det=SSE2(1,:).*(SSE2(3,:).*SSE2(6,:)-SSE2(5,:).^2) - ...
0087             SSE2(6,:).*SSE2(2,:).^2 + ...
0088             SSE2(4,:).*(SSE2(2,:).*SSE2(5,:)*2-SSE2(3,:).*SSE2(4,:));
0089         m1=SSE2(3,:).*SSE2(6,:)-SSE2(5,:).^2;
0090         m3=SSE2(1,:).*SSE2(6,:)-SSE2(4,:).^2;
0091         m6=SSE2(1,:).*SSE2(3,:)-SSE2(2,:).^2;
0092         m2=SSE2(4,:).*SSE2(5,:)-SSE2(2,:).*SSE2(6,:);
0093         m4=SSE2(2,:).*SSE2(5,:)-SSE2(3,:).*SSE2(4,:);
0094         m5=SSE2(2,:).*SSE2(4,:)-SSE2(1,:).*SSE2(5,:);
0095         a11=m1.*h(1,:)+m2.*h(2,:)+m4.*h(4,:);
0096         a12=m1.*h(2,:)+m2.*h(3,:)+m4.*h(5,:);
0097         a13=m1.*h(4,:)+m2.*h(5,:)+m4.*h(6,:);
0098         a21=m2.*h(1,:)+m3.*h(2,:)+m5.*h(4,:);
0099         a22=m2.*h(2,:)+m3.*h(3,:)+m5.*h(5,:);
0100         a23=m2.*h(4,:)+m3.*h(5,:)+m5.*h(6,:);
0101         a31=m4.*h(1,:)+m5.*h(2,:)+m6.*h(4,:);
0102         a32=m4.*h(2,:)+m5.*h(3,:)+m6.*h(5,:);
0103         a33=m4.*h(4,:)+m5.*h(5,:)+m6.*h(6,:);
0104         a0=-a11.*(a22.*a33-a23.*a32) ...
0105            +a12.*(a21.*a33-a23.*a31) ...
0106            -a13.*(a21.*a32-a22.*a31);
0107         a1=a22.*a33-a23.*a32 + ...
0108            a11.*a33-a13.*a31 + ...
0109            a11.*a22-a12.*a21;
0110         a2=-(a11+a22+a33);
0111         q=a1/3-a2.^2/9;
0112         r=(a1.*a2-3*a0)/6-a2.^3/27;
0113         s1=(r+sqrt(q.^3+r.^2)).^(1/3);
0114         z=zeros(3,v);
0115         z(1,:)=2*real(s1)-a2/3;
0116         z(2,:)=-real(s1)-a2/3+sqrt(3)*imag(s1);
0117         z(3,:)=-real(s1)-a2/3-sqrt(3)*imag(s1);
0118         z=sort(z,1,'descend');
0119         d=(df2/(df1-df2))./(det+(det<=0)).*(det>0);
0120         for j=1:l
0121             slm.t(j,:)=z(j,:).*d;
0122         end
0123     end
0124 end
0125 
0126 return
0127 end
0128