Attachment 'batch_spm5_stats_2ndlevel_mixed_ANOVA.m'
Download 1 % A quick hack for mixed (within+between subjects) ANOVAs
2 % see http://imaging.mrc-cbu.cam.ac.uk/imaging/SpmBatch5s. R Henson
3
4 spm_defaults
5 global defaults
6 global UFp
7 UFp = 0.001;
8 defaults.modality='FMRI';
9
10 % User input -----------------------------------------------------------------
11
12 ncon = 12; % number of conditions per subject
13 nbf = 1; % number of basis functions per condition
14 nsub = [12 15]; % number of subjects per group
15 cont = [34:45]; % contrast numbers from each subject's 1st-level
16
17 subjects_dir{1} = {
18 '/imaging/E20/intentionalParticipants/CBU080657/ExpRes'
19 '/imaging/E20/intentionalParticipants/CBU080658/ExpRes'
20 '/imaging/E20/intentionalParticipants/CBU080666/ExpRes'
21 '/imaging/E20/intentionalParticipants/CBU080673/ExpRes'
22 '/imaging/E20/intentionalParticipants/CBU080678/ExpRes'
23 '/imaging/E20/intentionalParticipants/CBU080694/ExpRes'
24 '/imaging/E20/intentionalParticipants/CBU080733/ExpRes'
25 '/imaging/E20/intentionalParticipants/CBU080735/ExpRes'
26 '/imaging/E20/intentionalParticipants/CBU080736/ExpRes'
27 '/imaging/E20/intentionalParticipants/CBU080737/ExpRes'
28 '/imaging/E20/intentionalParticipants/CBU080738/ExpRes'
29 '/imaging/E20/intentionalParticipants/CBU080780/ExpRes'
30 }
31 subjects_dir{2} = {
32 '/imaging/E20/incidentalParticipants/CBU080636/ExpRes'
33 '/imaging/E20/incidentalParticipants/CBU080659/ExpRes'
34 '/imaging/E20/incidentalParticipants/CBU080667/ExpRes'
35 '/imaging/E20/incidentalParticipants/CBU080671/ExpRes'
36 '/imaging/E20/incidentalParticipants/CBU080672/ExpRes'
37 '/imaging/E20/incidentalParticipants/CBU080697/ExpRes'
38 '/imaging/E20/incidentalParticipants/CBU080698/ExpRes'
39 '/imaging/E20/incidentalParticipants/CBU080699/ExpRes'
40 '/imaging/E20/incidentalParticipants/CBU080724/ExpRes'
41 '/imaging/E20/incidentalParticipants/CBU080726/ExpRes'
42 '/imaging/E20/incidentalParticipants/CBU080745/ExpRes'
43 '/imaging/E20/incidentalParticipants/CBU080747/ExpRes'
44 '/imaging/E20/incidentalParticipants/CBU080751/ExpRes'
45 '/imaging/E20/incidentalParticipants/CBU080778/ExpRes'
46 '/imaging/E20/incidentalParticipants/CBU080779/ExpRes'
47 };
48
49 group_dir = '/imaging/E20/CombIntInc/testphase/WeightedVisualStim';
50
51 pflag = 1; %whether you want figures of design matrix/nonsphericity before estimating
52
53 %-----------------------------------------------------------------
54 % Contrast name and numbers (for each condition/basis function)
55
56 totsub = sum(nsub); % (total number of subjects)
57 ngrp = length(nsub); % (number of groups)
58
59 ncol = ncon*nbf; % number of columns in X per group
60 nscan = sum(ncol.*nsub); % number of rows in X
61
62 try eval(sprintf('!mkdir %s',group_dir)); end
63 cd(group_dir);
64
65 % get image files names
66 P={}; cname={};
67 np=0; nc=0;
68 for g=1:ngrp
69 for n=1:ncol
70 for s=1:nsub(g)
71 np=np+1;
72 P{np} = fullfile(subjects_dir{g}{s},sprintf('con_%04d.img',cont(n)));
73 end
74 nc=nc+1;
75 cname{nc} = sprintf('grp %d, col %d',g,n);
76 end
77 end
78
79 for s=1:totsub % add subject effects
80 cname{end+1} = sprintf('subject %d',s);
81 end
82
83
84 %-Assemble SPM structure
85 %=======================================================================
86
87 SPM.nscan = nscan;
88 SPM.xY.P = P;
89 for i=1:SPM.nscan
90 SPM.xY.VY(i) = spm_vol(SPM.xY.P{i});
91 end
92
93 % Build design matrix (X), Indices (Ind) and NONSPHERICITY (vi) (inelegant, but gets there...!)
94
95 X=[]; Ind=[]; vi={};
96 nv=0; z=zeros(nscan,nscan); os=0;
97
98 for g=1:ngrp
99 ns = nsub(g);
100 nr = ncol*ns;
101 id = [1:ns]';
102
103 tmpX = [zeros(nr,ncol*(g-1)),...
104 kron(eye(ncol),ones(ns,1)),...
105 zeros(nr,ncol*(ngrp-g))];
106
107 % could add constants for group effects if wish
108
109 tmpX = [tmpX zeros(size(tmpX,1),sum(nsub(1:(g-1)))) kron(ones(ncol,1),eye(ns))];
110
111 if g>1
112 X = [X zeros(size(X,1),nsub(g)); tmpX];
113 else
114 X = tmpX;
115 end
116
117 % Indices for effects (unnecessary really)
118 Ind = [Ind; kron(ones(ncol,1),id),...
119 kron([1:ncol]',ones(ns,1)),...
120 ones(nr,1) ones(nr,1)*g];
121
122 % Nonsphericity
123
124 % unequal variances
125 for c1 = 1:ncol
126 nv = nv+1;
127 v = z;
128 v(os + (c1-1)*ns + id, os + (c1-1)*ns + id)=eye(ns);
129 vi{nv} = sparse(v);
130 end
131
132 % unequal covariances (within conditions; independent between groups)
133 for c1 = 1:(ncol-1)
134 for c2 = (c1+1):ncol
135 nv = nv+1;
136 v = z;
137 v( os + (c1-1)*ns + id, os + (c2-1)*ns + id )=eye(ns);
138 v( os + (c2-1)*ns + id, os + (c1-1)*ns + id )=eye(ns);
139 vi{nv} = sparse(v);
140 end
141 end
142 os = os + ncol*ns;
143 end
144
145 if pflag
146 figure,imagesc(X),colormap('gray')
147 figure,hold on,colormap('gray')
148 for pp=1:length(vi)
149 subplot(1,length(vi),pp)
150 imagesc(vi{pp})
151 end
152 pause
153 end
154
155 temp = [ones(nscan,1) kron((1:ncol)',ones(totsub,1)) kron(ones(ncol,1),(1:totsub)') ones(nscan,1)];
156 SPM.xX = struct(...
157 'X',X,...
158 'iH',[1:size(X,2)],'iC',zeros(1,0),'iB',zeros(1,0),'iG',zeros(1,0),...
159 'name',{cname},'I',temp,...
160 'sF',{{'repl' 'col' 'dummy' 'grp'}});
161 SPM.xC = [];
162 SPM.xGX = struct(...
163 'iGXcalc',1, 'sGXcalc','omit', 'rg',[],...
164 'iGMsca',9, 'sGMsca','<no grand Mean scaling>',...
165 'GM',0, 'gSF',ones(nscan,1),...
166 'iGC', 12, 'sGC', '(redundant: not doing AnCova)', 'gc',[],...
167 'iGloNorm',9, 'sGloNorm','<no global normalisation>');
168 SPM.xVi = struct('I',SPM.xX.I,'Vi',{vi} );
169 Mdes = struct(...
170 'Analysis_threshold', {'None (-Inf)'},...
171 'Implicit_masking', {'Yes: NaNs treated as missing'},...
172 'Explicit_masking', {'No'});
173 SPM.xM = struct(...
174 'T',-Inf,'TH',ones(nscan,1)*-Inf,...
175 'I',1,'VM',[],'xs',Mdes);
176 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)}};
177 SPM.xsDes = struct(...
178 'Design', {'1-way ANOVA'},...
179 'Global_calculation', {'omit'},...
180 'Grand_mean_scaling', {'<no grand Mean scaling>'},...
181 'Global_normalisation', {'<no global normalisation>'},...
182 'Parameters', Pdes);
183 save SPM SPM
184
185
186 % Estimate parameters
187 %===========================================================================
188 SPM = spm_spm(SPM);
189
190
191 % Effects of interest contrast?
192 %===========================================================================
193 %SPM = rmfield(SPM,'xCon'); cn=0;
194
195 cn = size(SPM.xCon,2);
196 c = detrend(eye(ncol*ngrp),0);
197 c = [c zeros(size(c,1),size(SPM.xX.X,2)-(ncol*ngrp))];
198 cname = 'Unwhitened effects of interest';
199 SPM.xCon(cn+1) = spm_FcUtil('Set',cname,'F','c',c',SPM.xX.xKXs);
200 spm_contrasts(SPM);
201
202 cn = size(SPM.xCon,2);
203 c = ones(1,ncol*ngrp);
204 c = [c zeros(size(c,1),size(SPM.xX.X,2)-(ncol*ngrp))];
205 cname = 'Activation vs Baseline';
206 SPM.xCon(cn+1) = spm_FcUtil('Set',cname,'T','c',c',SPM.xX.xKXs);
207 spm_contrasts(SPM);
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.You are not allowed to attach a file to this page.