1 % Batch for within-subject ANOVAs in SPM5       Rik Henson  Jan 06
   2 
   3 spm('defaults', 'FMRI')
   4 global defaults
   5 global UFp
   6 UFp=0.001;
   7 
   8 bwd = '/myhome';
   9 cd(bwd)
  10 
  11 subnam = [060001 060002 060003 060004 060005 060013 060014 060015 060016 060017 060018 060019 060023 060024 060025 060035 060036];
  12 nsub = length(subnam);
  13 
  14 ses = 'Ses1'; cons = [1:12]; ana = 'Mono12'; 
  15 
  16 ncon = length(cons);
  17 fnm = 'con';
  18 %fnm = 'beta';
  19 nsc = 1;	
  20 
  21 nscan = nsub*ncon;
  22 
  23 nwd = sprintf('%s/Group/%s/%s',bwd,ses,ana)
  24 eval(sprintf('!mkdir %s',nwd));
  25 cd(nwd);
  26 
  27 
  28 % get image files names
  29 P={};
  30 for c=1:ncon
  31 	con = cons(c);
  32 	for s=1:nsub
  33 	     	sub = subnam(s);
  34       	     	P{(c-1)*nsub+s} = sprintf('%s/CBU%06d/%s/Stats/%s_%04d.img',bwd,subnam(s),ses,fnm,con);
  35 	end 
  36 end
  37 
  38 %-Assemble SPM structure
  39 %=======================================================================
  40 clear SPM
  41 
  42 SPM.nscan = nscan;
  43 SPM.xY.P = P;
  44 
  45 for i=1:SPM.nscan
  46  SPM.xY.VY(i) = spm_vol(SPM.xY.P{i});
  47 end
  48 
  49 cname = cell(ncon,1);
  50 for c=1:ncon
  51 	cname{c} = sprintf('con %d',c);
  52 end
  53 for c=1:nsub
  54     cname{c+ncon} = sprintf('sub %d',c);
  55 end
  56 
  57 os = ones(nsub,1); oc = ones(ncon,1);
  58 SPM.xX = struct(...
  59         'X',[kron(eye(ncon),os) kron(oc,eye(nsub))] ,...
  60         'iH',[1:(ncon+nsub)],'iC',zeros(1,0),'iB',zeros(1,0),'iG',zeros(1,0),...
  61         'name',{cname},'I',[ones(nscan,1) kron([1:ncon]',os) kron(oc,[1:nsub]') ones(nscan,1)],...
  62         'sF',{{'repl'  'cond'  'subj'  ''}});
  63 
  64 SPM.xC = [];	
  65 
  66 SPM.xGX = struct(...
  67         'iGXcalc',1,    'sGXcalc','omit',                               'rg',[],...
  68         'iGMsca',9,     'sGMsca','<no grand Mean scaling>',...
  69         'GM',0,         'gSF',ones(nscan,1),...
  70         'iGC',  12,     'sGC',  '(redundant: not doing AnCova)',        'gc',[],...
  71         'iGloNorm',9,   'sGloNorm','<no global normalisation>');
  72 
  73 x=zeros(ncon);
  74 s=eye(nsub);
  75 nv=0;
  76 % first do the leading diagonal elements
  77 for i=1:ncon
  78         nv=nv+1;
  79         v=x;
  80         v(i,i)=1;
  81         vi{nv}=sparse(kron(v,s));
  82 end
  83 % then do the off-diagonals
  84 for i=1:(ncon-1)
  85         for j=(i+1):ncon
  86                 nv=nv+1;
  87                 v=x;
  88                 v(j,i)=1; v(i,j)=1;
  89                 vi{nv}=sparse(kron(v,s));
  90        end
  91 end
  92 
  93 if ncon==1 | nsc==0
  94  SPM.xVi = struct('iid',1, 'V',speye(nsub*ncon) );       
  95 else
  96  SPM.xVi = struct('iid',0, 'I',SPM.xX.I, 'Vi',{vi} );       
  97 end
  98 
  99 Mdes    = struct(...	
 100         'Analysis_threshold',   {'None (-Inf)'},...
 101         'Implicit_masking',     {'Yes: NaNs treated as missing'},...
 102         'Explicit_masking',     {'No'});
 103 
 104 SPM.xM  = struct(...
 105         'T',-Inf,'TH',ones(nscan,1)*-Inf,...
 106         'I',1,'VM',[],'xs',Mdes);
 107 
 108 Pdes    = {{sprintf('%d condition, +0 covariate, +%d block, +0 nuisance',ncon,nsub); sprintf('%d total, having %d degrees of freedom',ncon+nsub,ncon+nsub-1); sprintf('leaving %d degrees of freedom from %d images',nscan-ncon-nsub+1,nscan)}};
 109 
 110 Pdes{1}
 111 
 112 SPM.xsDes = struct(...
 113         'Design',               {'1-way ANOVA (within-subjects)'},...
 114         'Global_calculation',   {'omit'},...
 115         'Grand_mean_scaling',   {'<no grand Mean scaling>'},...
 116         'Global_normalisation', {'<no global normalisation>'},...
 117         'Parameters',           Pdes);
 118 
 119 save SPM SPM
 120 
 121 
 122 % Estimate parameters
 123 %===========================================================================
 124 SPM = spm_spm(SPM);
 125 
 126 
 127 SPM = rmfield(SPM,'xCon');
 128 
 129 c = eye(ncon);
 130 c = [c zeros(ncon,nsub)];
 131 cname       = 'Unwhitened EOI';
 132 SPM.xCon(1) = spm_FcUtil('Set',cname,'F','c',c',SPM.xX.xKXs);
 133 
 134 
 135 % F-contrasts 
 136 %-------------
 137 
 138 hemi = [1 1; 1 -1]; hemin = {'L+R';'L-R'};
 139 stim = [1 1; 1 -1]; stimn = {'F+H';'F-H'};
 140 reps = [1 1 1; 2 -1 -1]; repsn = {'1+2';'1-2'};
 141 
 142 cc=1
 143 for nh=1:size(hemi,1)
 144  for ns=1:size(stim,1)
 145   for nr=1:size(reps,1)
 146    c = kron(hemi(nh,:),stim(ns,:));
 147    c = kron(c,reps(nr,:));
 148    c = norm_cons(c);
 149    c = [c zeros(1,nsub)];
 150    cc=cc+1;
 151    cname = sprintf('%s X %s X %s',hemin{nh},stimn{ns},repsn{nr});
 152    SPM.xCon(cc) = spm_FcUtil('Set',cname,'F','c',c',SPM.xX.xKXs);
 153   end
 154  end
 155 end
 156 
 157  
 158 % T-contrasts 
 159 %-------------
 160 
 161 hemi = [1 1; 1 -1; -1 1]; hemin = {'L+R';'L-R';'R-L'};
 162 stim = [1 1; 1 -1; -1 1; 1 0; 0 1]; stimn = {'F+H';'F-H';'H-F';'F';'H'};
 163 reps = [1 1 1; 1 -1 -1; -1 1 1; 1 -1 0; -1 1 0];
 164 repsn = {'All';'1-2';'2-1';'1-S';'S-1'};
 165 
 166 %cc=1
 167 for nh=1:size(hemi,1)
 168  for ns=1:size(stim,1)
 169   for nr=1:size(reps,1)
 170    c = kron(hemi(nh,:),stim(ns,:));
 171    c = kron(c,reps(nr,:));
 172    c = norm_cons(c);
 173    c = [c zeros(1,nsub)];
 174    cc=cc+1;
 175    cname = sprintf('%s X %s X %s',hemin{nh},stimn{ns},repsn{nr});
 176    SPM.xCon(cc) = spm_FcUtil('Set',cname,'T','c',c',SPM.xX.xKXs);
 177   end
 178  end
 179 end
 180 
 181 spm_contrasts(SPM);