SurfStatLinMod

Fits linear mixed effects models to surface data and estimates resels.

 Usage: slm = SurfStatLinMod( Y, M [,surf [,niter [,thetalim [,drlim]]]]);

 Y        = n x v or n x v x k matrix of surface data, v=#vertices;
            n=#observations; k=#variates, or memory map of same.      
 M        = model formula of class term or random; or scalar or
            n x p design matrix of p regressors for the linear model.
 surf.tri = t x 3 matrix of triangle indices, 1-based, t=#triangles.
 or
 surf.lat = nx x ny x nz matrix, 1=in, 0=out, [nx,ny,nz]=size(volume). 
 The following arguments are not usually used:
 niter    = number of extra iterations of the Fisher scoring algorithm
            for fitting mixed effects models. Default 1.
 thetalim = lower limit on variance coefficients, in sd's. Default 0.01.
 drlim    = step of ratio of variance coefficients, in sd's. Default 0.1. 

 slm.X    = n x p design matrix.
 slm.V    = n x n x q variance matrix bases, normalised so that
            mean(diag(slm.V))=1. Absent if q=1 and slm.V=eye(n).
 slm.df   = degrees of freedom = n-rank(X).
 slm.coef = p x v x k matrix of coefficients of the linear model.
 slm.SSE  = k*(k+1)/2 x v matrix of sum of squares of errors
            whitened by slm.V*slm.r.
 slm.r    = (q-1) x v matrix of coefficients of the first (q-1)
            components of slm.V divided by their sum. 
            Coefficients are clamped to a minimum of 0.01 x sd.
 slm.dr   = (q-1) x nc vector of increments in slm.r = 0.1 x sd.
 If surf is not empty, returns:
 slm.resl = e x k matrix of sum over observations of squares of 
            differences of normalized residuals along each edge.
 slm.tri  = surf.tri,
 or
 slm.lat  = surf.lat.

0001 function slm = SurfStatLinMod( Y, M, surf, niter, thetalim, drlim );
0002 
0003 %Fits linear mixed effects models to surface data and estimates resels.
0004 %
0005 % Usage: slm = SurfStatLinMod( Y, M [,surf [,niter [,thetalim [,drlim]]]]);
0006 %
0007 % Y        = n x v or n x v x k matrix of surface data, v=#vertices;
0008 %            n=#observations; k=#variates, or memory map of same.
0009 % M        = model formula of class term or random; or scalar or
0010 %            n x p design matrix of p regressors for the linear model.
0011 % surf.tri = t x 3 matrix of triangle indices, 1-based, t=#triangles.
0012 % or
0013 % surf.lat = nx x ny x nz matrix, 1=in, 0=out, [nx,ny,nz]=size(volume).
0014 % The following arguments are not usually used:
0015 % niter    = number of extra iterations of the Fisher scoring algorithm
0016 %            for fitting mixed effects models. Default 1.
0017 % thetalim = lower limit on variance coefficients, in sd's. Default 0.01.
0018 % drlim    = step of ratio of variance coefficients, in sd's. Default 0.1.
0019 %
0020 % slm.X    = n x p design matrix.
0021 % slm.V    = n x n x q variance matrix bases, normalised so that
0022 %            mean(diag(slm.V))=1. Absent if q=1 and slm.V=eye(n).
0023 % slm.df   = degrees of freedom = n-rank(X).
0024 % slm.coef = p x v x k matrix of coefficients of the linear model.
0025 % slm.SSE  = k*(k+1)/2 x v matrix of sum of squares of errors
0026 %            whitened by slm.V*slm.r.
0027 % slm.r    = (q-1) x v matrix of coefficients of the first (q-1)
0028 %            components of slm.V divided by their sum.
0029 %            Coefficients are clamped to a minimum of 0.01 x sd.
0030 % slm.dr   = (q-1) x nc vector of increments in slm.r = 0.1 x sd.
0031 % If surf is not empty, returns:
0032 % slm.resl = e x k matrix of sum over observations of squares of
0033 %            differences of normalized residuals along each edge.
0034 % slm.tri  = surf.tri,
0035 % or
0036 % slm.lat  = surf.lat.
0037 
0038 maxchunk=2^20;
0039 if nargin<4 | isempty(niter)
0040     niter=1;
0041 end
0042 if nargin<5 | isempty(thetalim)
0043     thetalim=0.01;
0044 end
0045 if nargin<6 
0046     drlim=0.1;
0047 end
0048 
0049 if isnumeric(Y)
0050     s=size(Y);
0051     isnum=true;
0052 else
0053     Ym=Y;
0054     s=Ym.Format{2};
0055     if length(s)==2
0056         s=s([2 1]);
0057     else
0058         s=s([3 1 2]);
0059     end
0060     [PATHSTR,NAME,EXT]=fileparts(Ym.Filename);
0061     filenameResid=fullfile(PATHSTR,['Resid.mmm']);
0062     copyfile(Ym.Filename,filenameResid);
0063     YmResid=memmapfile(filenameResid,'Format',{'single' s 'Resid'},'Writable',true);
0064     isnum=false;
0065 end
0066 n=s(1);
0067 v=s(2);
0068 if length(s)==2
0069     k=1;
0070 else
0071     k=s(3);
0072 end
0073 
0074 if isa(M,'random')
0075     [slm.X,V]=double(M);
0076     [n2 q]=size(V);
0077     II=reshape(eye(n),n^2,1);
0078     r=II-V*(pinv(V)*II);
0079     if mean(r.^2)>eps
0080         warning('Did you forget an error term, I? :-)');
0081     end
0082     if q>1 | ((q==1) & sum(abs(II-V))>0)
0083         slm.V=reshape(V,[n n q]);
0084     end
0085     clear V II r;
0086 else
0087     q=1;
0088     if isa(M,'term')
0089         slm.X=double(M);
0090     else
0091         if prod(size(M))>1
0092             warning('If you don''t convert vectors to terms you can get unexpected results :-(')
0093         end
0094         slm.X=M;
0095     end
0096     if size(slm.X,1)==1
0097         slm.X=repmat(slm.X,n,1);
0098     end
0099 end
0100 
0101 pinvX=pinv(slm.X);
0102 r=ones(n,1)-slm.X*(pinvX*ones(n,1));
0103 if mean(r.^2)>eps
0104     warning('Did you forget a constant term? :-)');
0105 end
0106 
0107 p=size(slm.X,2);
0108 slm.df=n-rank(slm.X);
0109 slm.coef=zeros(p,v,k);
0110 k2=k*(k+1)/2;
0111 slm.SSE=zeros(k2,v);
0112 
0113 if isnum
0114     nc=1;
0115     chunk=v;
0116 else
0117     nc=ceil(v*n*k/maxchunk);
0118     chunk=ceil(v/nc);
0119 end
0120 if ~isnum
0121     fprintf(1,'%s',[num2str(round(v*n*k*4/2^20)) ' Mb to fit, % remaining: 100 ']);
0122     n10=floor(n/10);
0123 end
0124 for ic=1:nc
0125     if ~isnum & rem(ic,n10)==0 
0126         fprintf(1,'%s',[num2str(round(100*(1-ic/nc))) ' ']);
0127     end
0128     v1=1+(ic-1)*chunk;
0129     v2=min(v1+chunk-1,v);
0130     vc=v2-v1+1;
0131     if isnum       
0132         Y=double(Y);
0133     else
0134         if length(s)==2
0135             Y=double(Ym.Data(1).Data(v1:v2,:)');
0136         else
0137             Y=double(permute(Ym.Data(1).Data(v1:v2,:,:),[3 1 2]));
0138         end
0139     end
0140     if k==1
0141         if q==1
0142             %% fixed effects
0143             if ~isfield(slm,'V')
0144                 coef=pinvX*Y;
0145                 Y=Y-slm.X*coef;
0146             else
0147                 if ic==1
0148                     slm.V=slm.V/mean(diag(slm.V));
0149                     Vmh=inv(chol(slm.V)');
0150                 end
0151                 coef=(pinv(Vmh*slm.X)*Vmh)*Y;
0152                 Y=Vmh*Y-(Vmh*slm.X)*coef;
0153             end
0154             SSE=sum(Y.^2);
0155         else
0156             %% mixed effects
0157             if ic==1
0158                 q1=q-1;
0159                 for j=1:q
0160                     slm.V(:,:,j)=slm.V(:,:,j)/mean(diag(slm.V(:,:,j)));;
0161                 end
0162                 slm.r=zeros(q1,v);
0163                 slm.dr=zeros(q1,nc);
0164             end
0165             coef=zeros(p,vc);
0166             SSE=zeros(1,vc);
0167 
0168             %% start Fisher scoring algorithm
0169             E=zeros(q,vc);
0170             RVV=zeros([n n q]);
0171             R=eye(n)-slm.X*pinv(slm.X);
0172             for j=1:q
0173                 RV=R*slm.V(:,:,j);
0174                 E(j,:)=sum(Y.*((RV*R)*Y));
0175                 RVV(:,:,j)=RV;
0176             end
0177             M=zeros(q);
0178             for j1=1:q
0179                 for j2=j1:q
0180                     M(j1,j2)=sum(sum(RVV(:,:,j1).*(RVV(:,:,j2)')));
0181                     M(j2,j1)=M(j1,j2);
0182                 end
0183             end
0184             theta=pinv(M)*E;
0185 
0186             tlim=sqrt(2*diag(pinv(M)))*sum(theta)*thetalim;
0187             theta=theta.*(theta>=tlim)+tlim.*(theta<tlim);
0188             r=theta(1:q1,:)./repmat(sum(theta),q1,1);
0189 
0190             Vt=2*pinv(M);
0191             m1=diag(Vt);
0192             m2=2*sum(Vt)';
0193             Vr=m1(1:q1)-m2(1:q1).*mean(slm.r,2)+sum(Vt(:))*mean(r.^2,2);
0194             dr=sqrt(Vr)*drlim;
0195 
0196             %% Exrtra Fisher scoring iterations
0197             for iter=1:niter
0198                 irs=round(r.*repmat(1./dr,1,vc));
0199                 [ur,ir,jr]=unique(irs','rows');
0200                 nr=size(ur,1);
0201                 for ir=1:nr
0202                     iv=(jr==ir);
0203                     rv=mean(r(:,iv),2);
0204                     V=(1-sum(rv))*slm.V(:,:,q);
0205                     for j=1:q1
0206                         V=V+rv(j)*slm.V(:,:,j);
0207                     end
0208                     Vinv=inv(V);
0209                     VinvX=Vinv*slm.X;
0210                     G=pinv(slm.X'*VinvX)*(VinvX');
0211                     R=Vinv-VinvX*G;
0212                     E=zeros(q,sum(iv));
0213                     for j=1:q
0214                         RV=R*slm.V(:,:,j);
0215                         E(j,:)=sum(Y(:,iv).*((RV*R)*Y(:,iv)));
0216                         RVV(:,:,j)=RV;
0217                     end
0218                     for j1=1:q
0219                         for j2=j1:q
0220                             M(j1,j2)=sum(sum(RVV(:,:,j1).*(RVV(:,:,j2)')));
0221                             M(j2,j1)=M(j1,j2);
0222                         end
0223                     end
0224                     thetav=pinv(M)*E;
0225                     tlim=sqrt(2*diag(pinv(M)))*sum(thetav)*thetalim;
0226                     theta(:,iv)=thetav.*(thetav>=tlim)+tlim.*(thetav<tlim);
0227                 end
0228                 r=theta(1:q1,:)./(ones(q1,1)*sum(theta));
0229             end
0230 
0231             %% finish Fisher scoring
0232             irs=round(r.*repmat(1./dr,1,vc));
0233             [ur,ir,jr]=unique(irs','rows');
0234             nr=size(ur,1);
0235             for ir=1:nr
0236                 iv=(jr==ir);
0237                 rv=mean(r(:,iv),2);
0238                 V=(1-sum(rv))*slm.V(:,:,q);
0239                 for j=1:q1
0240                     V=V+rv(j)*slm.V(:,:,j);
0241                 end
0242                 Vmh=inv(chol(V)');
0243                 VmhX=Vmh*slm.X;
0244                 G=pinv(VmhX'*VmhX)*(VmhX')*Vmh;
0245                 coef(:,iv)=G*Y(:,iv);
0246                 R=Vmh-VmhX*G;
0247                 Y(:,iv)=R*Y(:,iv);
0248                 SSE(iv)=sum(Y(:,iv).^2);
0249             end
0250             slm.r(:,v1:v2)=r;
0251             slm.dr(:,ic)=dr;
0252         end
0253     else
0254         %% multivariate
0255         if q>1
0256             error('Multivariate mixed effects models not yet implemented :-(');
0257             return
0258         end
0259 
0260         if ~isfield(slm,'V')
0261             X=slm.X;
0262         else
0263             if ic==1
0264                 slm.V=slm.V/mean(diag(slm.V));
0265                 Vmh=inv(chol(slm.V))';
0266                 X=Vmh*slm.X;
0267                 pinvX=pinv(X);
0268             end
0269             for j=1:k
0270                 Y(:,:,j)=Vmh*Y(:,:,j);
0271             end
0272         end
0273 
0274         coef=zeros(p,vc,k);
0275         for j=1:k
0276             coef(:,:,j)=pinvX*Y(:,:,j);
0277             Y(:,:,j)=Y(:,:,j)-X*coef(:,:,j);
0278         end
0279         k2=k*(k+1)/2;
0280         SSE=zeros(k2,vc);
0281         j=0;
0282         for j1=1:k
0283             for j2=1:j1
0284                 j=j+1;
0285                 SSE(j,:)=sum(Y(:,:,j1).*Y(:,:,j2));
0286             end
0287         end
0288     end
0289     slm.coef(:,v1:v2,:)=coef;
0290     slm.SSE(:,v1:v2)=SSE;
0291     if ~isnum
0292         YmResid.Data(1).Resid(:,v1:v2,:)=single(Y);
0293     end
0294 end
0295 if ~isnum
0296     fprintf(1,'%s\n','Done');
0297 end
0298 
0299 %% resels
0300 if nargin<3 | isempty(surf)
0301     return
0302 end
0303 
0304 if isfield(surf,'tri')
0305     slm.tri=surf.tri;
0306 end
0307 if isfield(surf,'lat')
0308     slm.lat=surf.lat;
0309 end
0310 
0311 edg=SurfStatEdg(surf);
0312 
0313 e=size(edg,1);
0314 e1=edg(:,1);
0315 e2=edg(:,2);
0316 clear edg
0317 slm.resl=zeros(e,k);
0318 if ~isnum
0319     fprintf(1,'%s',[num2str(n) ' x ' num2str(k) ' surfaces to resel, % remaining: 100 ']);
0320     n10=floor(n*k/10);
0321 end
0322 for j=1:k
0323     jj=j*(j+1)/2;
0324     normr=sqrt(slm.SSE(jj,:));
0325     s=0;
0326     for i=1:n
0327         if ~isnum & rem((j-1)*n+i,n10)==0
0328             fprintf(1,'%s',[num2str(round(100*(1-((j-1)*n+i)/(n*k)))) ' ']);
0329         end
0330         if isnum
0331             u=Y(i,:,j)./normr;
0332         else
0333             u=double(YmResid.Data(1).Resid(i,:,j))./normr;
0334         end
0335         s=s+(u(e1)-u(e2)).^2;
0336     end
0337     slm.resl(:,j)=s;
0338 end
0339 if ~isnum
0340     clear YmResid
0341     delete(filenameResid);
0342     fprintf(1,'%s\n','Done');
0343 end
0344 
0345 return
0346 end