% A quick hack for mixed (within+between subjects) ANOVAs % see http://imaging.mrc-cbu.cam.ac.uk/imaging/SpmBatch5s. R Henson spm_defaults global defaults global UFp UFp = 0.001; defaults.modality='FMRI'; % User input ----------------------------------------------------------------- ncon = 12; % number of conditions per subject nbf = 1; % number of basis functions per condition nsub = [12 15]; % number of subjects per group cont = [34:45]; % contrast numbers from each subject's 1st-level subjects_dir{1} = { '/imaging/E20/intentionalParticipants/CBU080657/ExpRes' '/imaging/E20/intentionalParticipants/CBU080658/ExpRes' '/imaging/E20/intentionalParticipants/CBU080666/ExpRes' '/imaging/E20/intentionalParticipants/CBU080673/ExpRes' '/imaging/E20/intentionalParticipants/CBU080678/ExpRes' '/imaging/E20/intentionalParticipants/CBU080694/ExpRes' '/imaging/E20/intentionalParticipants/CBU080733/ExpRes' '/imaging/E20/intentionalParticipants/CBU080735/ExpRes' '/imaging/E20/intentionalParticipants/CBU080736/ExpRes' '/imaging/E20/intentionalParticipants/CBU080737/ExpRes' '/imaging/E20/intentionalParticipants/CBU080738/ExpRes' '/imaging/E20/intentionalParticipants/CBU080780/ExpRes' } subjects_dir{2} = { '/imaging/E20/incidentalParticipants/CBU080636/ExpRes' '/imaging/E20/incidentalParticipants/CBU080659/ExpRes' '/imaging/E20/incidentalParticipants/CBU080667/ExpRes' '/imaging/E20/incidentalParticipants/CBU080671/ExpRes' '/imaging/E20/incidentalParticipants/CBU080672/ExpRes' '/imaging/E20/incidentalParticipants/CBU080697/ExpRes' '/imaging/E20/incidentalParticipants/CBU080698/ExpRes' '/imaging/E20/incidentalParticipants/CBU080699/ExpRes' '/imaging/E20/incidentalParticipants/CBU080724/ExpRes' '/imaging/E20/incidentalParticipants/CBU080726/ExpRes' '/imaging/E20/incidentalParticipants/CBU080745/ExpRes' '/imaging/E20/incidentalParticipants/CBU080747/ExpRes' '/imaging/E20/incidentalParticipants/CBU080751/ExpRes' '/imaging/E20/incidentalParticipants/CBU080778/ExpRes' '/imaging/E20/incidentalParticipants/CBU080779/ExpRes' }; group_dir = '/imaging/E20/CombIntInc/testphase/WeightedVisualStim'; pflag = 1; %whether you want figures of design matrix/nonsphericity before estimating %----------------------------------------------------------------- % Contrast name and numbers (for each condition/basis function) totsub = sum(nsub); % (total number of subjects) ngrp = length(nsub); % (number of groups) ncol = ncon*nbf; % number of columns in X per group nscan = sum(ncol.*nsub); % number of rows in X try eval(sprintf('!mkdir %s',group_dir)); end cd(group_dir); % get image files names P={}; cname={}; np=0; nc=0; for g=1:ngrp for n=1:ncol for s=1:nsub(g) np=np+1; P{np} = fullfile(subjects_dir{g}{s},sprintf('con_%04d.img',cont(n))); end nc=nc+1; cname{nc} = sprintf('grp %d, col %d',g,n); end end for s=1:totsub % add subject effects cname{end+1} = sprintf('subject %d',s); end %-Assemble SPM structure %======================================================================= SPM.nscan = nscan; SPM.xY.P = P; for i=1:SPM.nscan SPM.xY.VY(i) = spm_vol(SPM.xY.P{i}); end % Build design matrix (X), Indices (Ind) and NONSPHERICITY (vi) (inelegant, but gets there...!) X=[]; Ind=[]; vi={}; nv=0; z=zeros(nscan,nscan); os=0; for g=1:ngrp ns = nsub(g); nr = ncol*ns; id = [1:ns]'; tmpX = [zeros(nr,ncol*(g-1)),... kron(eye(ncol),ones(ns,1)),... zeros(nr,ncol*(ngrp-g))]; % could add constants for group effects if wish tmpX = [tmpX zeros(size(tmpX,1),sum(nsub(1:(g-1)))) kron(ones(ncol,1),eye(ns))]; if g>1 X = [X zeros(size(X,1),nsub(g)); tmpX]; else X = tmpX; end % Indices for effects (unnecessary really) Ind = [Ind; kron(ones(ncol,1),id),... kron([1:ncol]',ones(ns,1)),... ones(nr,1) ones(nr,1)*g]; % Nonsphericity % unequal variances for c1 = 1:ncol nv = nv+1; v = z; v(os + (c1-1)*ns + id, os + (c1-1)*ns + id)=eye(ns); vi{nv} = sparse(v); end % unequal covariances (within conditions; independent between groups) for c1 = 1:(ncol-1) for c2 = (c1+1):ncol nv = nv+1; v = z; v( os + (c1-1)*ns + id, os + (c2-1)*ns + id )=eye(ns); v( os + (c2-1)*ns + id, os + (c1-1)*ns + id )=eye(ns); vi{nv} = sparse(v); end end os = os + ncol*ns; end if pflag figure,imagesc(X),colormap('gray') figure,hold on,colormap('gray') for pp=1:length(vi) subplot(1,length(vi),pp) imagesc(vi{pp}) end pause end temp = [ones(nscan,1) kron((1:ncol)',ones(totsub,1)) kron(ones(ncol,1),(1:totsub)') ones(nscan,1)]; SPM.xX = struct(... 'X',X,... 'iH',[1:size(X,2)],'iC',zeros(1,0),'iB',zeros(1,0),'iG',zeros(1,0),... 'name',{cname},'I',temp,... 'sF',{{'repl' 'col' 'dummy' 'grp'}}); SPM.xC = []; SPM.xGX = struct(... 'iGXcalc',1, 'sGXcalc','omit', 'rg',[],... 'iGMsca',9, 'sGMsca','',... 'GM',0, 'gSF',ones(nscan,1),... 'iGC', 12, 'sGC', '(redundant: not doing AnCova)', 'gc',[],... 'iGloNorm',9, 'sGloNorm',''); SPM.xVi = struct('I',SPM.xX.I,'Vi',{vi} ); Mdes = struct(... 'Analysis_threshold', {'None (-Inf)'},... 'Implicit_masking', {'Yes: NaNs treated as missing'},... 'Explicit_masking', {'No'}); SPM.xM = struct(... 'T',-Inf,'TH',ones(nscan,1)*-Inf,... 'I',1,'VM',[],'xs',Mdes); Pdes = {{sprintf('%d condition, +0 covariate, +0 block, +0 nuisance',ncon); sprintf('%d total, having %d degrees of freedom',ncon,ncon); sprintf('leaving %d degrees of freedom from %d images',nscan-ncon,nscan)}}; SPM.xsDes = struct(... 'Design', {'1-way ANOVA'},... 'Global_calculation', {'omit'},... 'Grand_mean_scaling', {''},... 'Global_normalisation', {''},... 'Parameters', Pdes); save SPM SPM % Estimate parameters %=========================================================================== SPM = spm_spm(SPM); % Effects of interest contrast? %=========================================================================== %SPM = rmfield(SPM,'xCon'); cn=0; cn = size(SPM.xCon,2); c = detrend(eye(ncol*ngrp),0); c = [c zeros(size(c,1),size(SPM.xX.X,2)-(ncol*ngrp))]; cname = 'Unwhitened effects of interest'; SPM.xCon(cn+1) = spm_FcUtil('Set',cname,'F','c',c',SPM.xX.xKXs); spm_contrasts(SPM); cn = size(SPM.xCon,2); c = ones(1,ncol*ngrp); c = [c zeros(size(c,1),size(SPM.xX.X,2)-(ncol*ngrp))]; cname = 'Activation vs Baseline'; SPM.xCon(cn+1) = spm_FcUtil('Set',cname,'T','c',c',SPM.xX.xKXs); spm_contrasts(SPM);