T statistics for a contrast in a univariate or multivariate model.
Usage: slm = SurfStatT( slm, contrast );
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. If absent, assumes 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.
slm.r = (q-1) x v matrix of coefficients of the first (q-1)
components of slm.V divided by their sum.
slm.dr = (q-1) x 1 vector of increments in slm.r.
contrast = n x 1 vector of contrasts in the observations, i.e.
= slm.X*slm.c', where slm.c is a contrast in slm.coef, or
= 1 x p vector of slm.c, padded with 0's if length(contrast)<p.
slm.c = 1 x p vector of contrasts in coefficents of the linear model.
slm.k = k=#variates.
slm.ef = k x v matrix of effects.
slm.sd = k x v matrix of standard deviations of the effects.
slm.t = 1 x v vector of T = ef/sd if k=1, or Hotelling's T if k=2 or 3,
defined as the maximum T over all linear combinations of the k
variates; k>3 is not programmed yet.
slm.dfs = 1 x v vector of effective degrees of freedom. Absent if q=1.
Note that the contrast in the observations is used to determine the
intended contrast in the model coefficients, slm.c. However there is some
ambiguity in this when the model contains redundant terms. An example of
such a model is 1 + Gender (Gender by itself does not contain redundant
terms). Only one of the ambiguous contrasts is estimable (i.e. has slm.sd
< Inf), and this is the one chosen, though it may not be the contrast
that you intended. To check this, compare the contrast in the
coefficients slm.c to the actual design matrix in slm.X. Note that the
redundant columns of the design matrix have weights given by the rows of
null(slm.X,'r')'0001 function slm = SurfStatT( slm, contrast ); 0002 0003 %T statistics for a contrast in a univariate or multivariate model. 0004 % 0005 % Usage: slm = SurfStatT( slm, contrast ); 0006 % 0007 % slm.X = n x p design matrix. 0008 % slm.V = n x n x q variance matrix bases, normalised so that 0009 % mean(diag(slm.V))=1. If absent, assumes q=1 and slm.V=eye(n). 0010 % slm.df = degrees of freedom = n-rank(X). 0011 % slm.coef = p x v x k matrix of coefficients of the linear model. 0012 % slm.SSE = k*(k+1)/2 x v matrix of sum of squares of errors. 0013 % slm.r = (q-1) x v matrix of coefficients of the first (q-1) 0014 % components of slm.V divided by their sum. 0015 % slm.dr = (q-1) x 1 vector of increments in slm.r. 0016 % contrast = n x 1 vector of contrasts in the observations, i.e. 0017 % = slm.X*slm.c', where slm.c is a contrast in slm.coef, or 0018 % = 1 x p vector of slm.c, padded with 0's if length(contrast)<p. 0019 % 0020 % slm.c = 1 x p vector of contrasts in coefficents of the linear model. 0021 % slm.k = k=#variates. 0022 % slm.ef = k x v matrix of effects. 0023 % slm.sd = k x v matrix of standard deviations of the effects. 0024 % slm.t = 1 x v vector of T = ef/sd if k=1, or Hotelling's T if k=2 or 3, 0025 % defined as the maximum T over all linear combinations of the k 0026 % variates; k>3 is not programmed yet. 0027 % slm.dfs = 1 x v vector of effective degrees of freedom. Absent if q=1. 0028 % 0029 % Note that the contrast in the observations is used to determine the 0030 % intended contrast in the model coefficients, slm.c. However there is some 0031 % ambiguity in this when the model contains redundant terms. An example of 0032 % such a model is 1 + Gender (Gender by itself does not contain redundant 0033 % terms). Only one of the ambiguous contrasts is estimable (i.e. has slm.sd 0034 % < Inf), and this is the one chosen, though it may not be the contrast 0035 % that you intended. To check this, compare the contrast in the 0036 % coefficients slm.c to the actual design matrix in slm.X. Note that the 0037 % redundant columns of the design matrix have weights given by the rows of 0038 % null(slm.X,'r')' 0039 0040 [n,p]=size(slm.X); 0041 contrast=double(contrast); 0042 pinvX=pinv(slm.X); 0043 if length(contrast)<=p 0044 c=[contrast zeros(1,p-length(contrast))]'; 0045 if sum((null(slm.X)'*c).^2)/sum(c.^2)>eps 0046 error('Contrast is not estimable :-('); 0047 return 0048 end 0049 else 0050 c=pinvX*contrast; 0051 r=contrast-slm.X*c; 0052 if sum(r(:).^2)/sum(contrast(:).^2)>eps 0053 warning('Contrast is not in the model :-('); 0054 end 0055 end 0056 slm.c=c'; 0057 0058 slm.df=slm.df(length(slm.df)); 0059 if ndims(slm.coef)==2 0060 slm.k=1; 0061 if ~isfield(slm,'r') 0062 %% fixed effect 0063 if isfield(slm,'V') 0064 Vmh=inv(chol(slm.V)'); 0065 pinvX=pinv(Vmh*slm.X); 0066 end 0067 Vc=sum((c'*pinvX).^2,2); 0068 else 0069 %% mixed effect 0070 [q1,v]=size(slm.r); 0071 q=q1+1; 0072 0073 nc=size(slm.dr,2); 0074 chunk=ceil(v/nc); 0075 irs=zeros(q1,v); 0076 for ic=1:nc 0077 v1=1+(ic-1)*chunk; 0078 v2=min(v1+chunk-1,v); 0079 vc=v2-v1+1; 0080 irs(:,v1:v2)=round(slm.r(:,v1:v2).*repmat(1./slm.dr(:,ic),1,vc)); 0081 end 0082 [ur,ir,jr]=unique(irs','rows'); 0083 nr=size(ur,1); 0084 slm.dfs=zeros(1,v); 0085 Vc=zeros(1,v); 0086 for ir=1:nr 0087 iv=(jr==ir); 0088 rv=mean(slm.r(:,iv),2); 0089 V=(1-sum(rv))*slm.V(:,:,q); 0090 for j=1:q1 0091 V=V+rv(j)*slm.V(:,:,j); 0092 end 0093 Vinv=inv(V); 0094 VinvX=Vinv*slm.X; 0095 Vbeta=pinv(slm.X'*VinvX); 0096 G=Vbeta*(VinvX'); 0097 Gc=G'*c; 0098 R=Vinv-VinvX*G; 0099 E=zeros(q,1); 0100 for j=1:q 0101 E(j)=Gc'*slm.V(:,:,j)*Gc; 0102 RVV(:,:,j)=R*slm.V(:,:,j); 0103 end 0104 for j1=1:q 0105 for j2=j1:q 0106 M(j1,j2)=sum(sum(RVV(:,:,j1).*(RVV(:,:,j2)'))); 0107 M(j2,j1)=M(j1,j2); 0108 end 0109 end 0110 vc=c'*Vbeta*c; 0111 iv=(ir==jr); 0112 Vc(iv)=vc; 0113 slm.dfs(iv)=vc^2/(E'*pinv(M)*E); 0114 end 0115 end 0116 slm.ef=c'*slm.coef; 0117 slm.sd=sqrt(Vc.*slm.SSE/slm.df); 0118 slm.t=slm.ef./(slm.sd+(slm.sd<=0)).*(slm.sd>0); 0119 else 0120 %% multivariate 0121 [p,v,k]=size(slm.coef); 0122 slm.k=k; 0123 slm.ef=zeros(k,v); 0124 for j=1:k 0125 slm.ef(j,:)=c'*slm.coef(:,:,j); 0126 end 0127 j=1:k; 0128 jj=j.*(j+1)/2; 0129 vf=sum((c'*pinvX).^2,2)/slm.df; 0130 slm.sd=sqrt(vf*slm.SSE(jj,:)); 0131 if k==2 0132 det=slm.SSE(1,:).*slm.SSE(3,:)-slm.SSE(2,:).^2; 0133 slm.t=slm.ef(1,:).^2.*slm.SSE(3,:) + ... 0134 slm.ef(2,:).^2.*slm.SSE(1,:) - ... 0135 2*slm.ef(1,:).*slm.ef(2,:).*slm.SSE(2,:); 0136 end 0137 if k==3 0138 det=slm.SSE(1,:).*(slm.SSE(3,:).*slm.SSE(6,:)-slm.SSE(5,:).^2) - ... 0139 slm.SSE(6,:).*slm.SSE(2,:).^2 + ... 0140 slm.SSE(4,:).*(slm.SSE(2,:).*slm.SSE(5,:)*2-slm.SSE(3,:).*slm.SSE(4,:)); 0141 slm.t= slm.ef(1,:).^2.*(slm.SSE(3,:).*slm.SSE(6,:)-slm.SSE(5,:).^2); 0142 slm.t=slm.t+slm.ef(2,:).^2.*(slm.SSE(1,:).*slm.SSE(6,:)-slm.SSE(4,:).^2); 0143 slm.t=slm.t+slm.ef(3,:).^2.*(slm.SSE(1,:).*slm.SSE(3,:)-slm.SSE(2,:).^2); 0144 slm.t=slm.t+2*slm.ef(1,:).*slm.ef(2,:).*(slm.SSE(4,:).*slm.SSE(5,:)-slm.SSE(2,:).*slm.SSE(6,:)); 0145 slm.t=slm.t+2*slm.ef(1,:).*slm.ef(3,:).*(slm.SSE(2,:).*slm.SSE(5,:)-slm.SSE(3,:).*slm.SSE(4,:)); 0146 slm.t=slm.t+2*slm.ef(2,:).*slm.ef(3,:).*(slm.SSE(2,:).*slm.SSE(4,:)-slm.SSE(1,:).*slm.SSE(5,:)); 0147 end 0148 if k>3 0149 warning('Hotelling''s T for k>3 not programmed yet'); 0150 return 0151 end 0152 slm.t=slm.t./(det+(det<=0)).*(det>0)/vf; 0153 slm.t=sqrt(slm.t+(slm.t<=0)).*(slm.t>0); 0154 end 0155 0156 return 0157 end 0158 0159 0160