-
Notifications
You must be signed in to change notification settings - Fork 10
/
snpm_pi_Corr1S.m
284 lines (262 loc) · 10.5 KB
/
snpm_pi_Corr1S.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
% Mfile snpm_pi_Corr1S
% SnPM PlugIn design module - SingleSubj simple correlation
% SingleSub: Simple Regression (correlation); single covariate of interest
% FORMAT snpm_pi_Corr1S
%
% See body of snpm_ui for definition of PlugIU interface.
%_______________________________________________________________________
%
% snpm_pi_Corr1S is a PlugIn for the SnPM design set-up program,
% creating design and permutation matrix appropriate for single-
% subject, correlation design.
%
%-Number of permutations
%=======================================================================
% There are nScan! (nScan factoral) possible permutations, where nScan
% is the number of scans of the subject. You can compute this using
% the gamma function in Matlab: nScan! is gamma(nScan+1); or by direct
% computation as prod(1:nScan)
%
%-Prompts
%=======================================================================
%
% 'Select scans in time order': Enter the scans to be analyzed.
% It is important to input the scans in time order so that temporal
% effects can be accounted for.
%
% 'Enter covariate values': These are the values of the experimental
% covariate.
%
% 'Size of exchangability block': This is the number of adjacent
% scans that you believe to be exchangeable. The most common cause
% of nonexchangeability is a temporal confound (e.g. while 12 adjacent
% scans might not be free from a temporal effect, 4 adjacent scans
% could be regarded as having neglible temporal effect). See snpm.man
% for more information and references.
%
% '### Perms. Use approx. test?': If there are a large number of
% permutations it may not be necessary (or possible!) to compute
% all the permutations. A common guideline is that 10,000 permutations
% are sufficient to characterize the permutation distribution well.
% More permutations will probably not change the results much.
% If you answer yes you will be prompted for the number...
% '# perms. to use?'
%
%
%-Variable "decoder" - This PlugIn supplies the following:
%=======================================================================
% - core -
% P - string matrix of Filenames corresponding to observations
% iGloNorm - Global normalisation code, or allowable codes
% - Names of columns of design matrix subpartitions
% PiCond - Permuted conditions matrix, one labelling per row, actual
% labelling on first row
% sPiCond - String describing permutations in PiCond
% sHCform - String for computation of HC design matrix partitions
% permutations indexed by perm in snpm_cp
% CONT - single contrast for examination, a row vector
% sDesign - String defining the design
% sDesSave - String of PlugIn variables to save to cfg file
%
% - design -
% C,Cnames - Condition partition of design matrix, & effect names
% B,Bnames - Block partition (constant term), & effect names
%
% - extra -
% CovInt - Specified covariate of interest
% Xblk - Size of exchangability block
%_______________________________________________________________________
% Copyright (C) 2013 The University of Warwick
% Id: snpm_pi_Corr1S.m SnPM13 2013/10/12
% Thomas Nichols & Andrew Holmes
%-----------------------------functions-called------------------------
% spm_DesMtx
% spm_select
% spm_input
%-----------------------------functions-called------------------------
%-Initialisation
%-----------------------------------------------------------------------
iGloNorm = '123'; % Allowable Global norm. codes
sDesSave = 'CovInt Xblk'; % PlugIn variables to save in cfg file
if snpm_get_defaults('shuffle_seed')
% Shuffle seed of random number generator
try
rng('shuffle');
catch
% Old syntax
rand('seed',sum(100*clock));
end
end
%-Get filenames of scans
%-----------------------------------------------------------------------
P = strvcat(job.P);%spm_select(Inf,'image','Select scans in time order');
nScan = size(P,1);
%-Get covariate
%-----------------------------------------------------------------------
CovInt = [];
while ~all(size(CovInt)==[nScan,1])
CovInt = job.CovInt(:);% spm_input(sprintf('Enter covariate values [%d]',nScan),'+1');
%CovInt = CovInt(:);
end
%-Centre covariate if required
%-----------------------------------------------------------------------
%**** bCCovInt = (spm_input('Centre this covariate ?','+0','y/n')=='y');
bCCovInt = 1; % Always centre covariate
%-Get confounding covariates
%-----------------------------------------------------------------------
G = []; Gnames = ''; Gc = []; Gcnames = ''; q = nScan;
if numel(job.cov) > 0 %isfield(job.covariate,'cov_Val')
for i = 1:numel(job.cov)
d = job.cov(i).c;
if (size(d,1) == 1)
d = d';
end
nGcs = size(Gc,2);
if size(d,1) ~= q
error('SnPM:InvalidCovariate', sprintf('Covariate [%d,1] does not match number of scans [%d]',...
size(job.cov(i).c,1),nScan))
else
%-Save raw covariates for printing later on
Gc = [Gc,d];
% Center
d = d - ones(q,1)*mean(d); str='';
G = [G, d];
dnames = job.cov(i).cname;
% dnames = [str,'ConfCov#',int2str(nGcs+1)];
% for j = nGcs+1:nGcs+size(d,1)
% dnames = str2mat(dnames,['ConfCov#',int2str(j)]);
% end
Gcnames = str2mat(Gcnames,dnames);
end
end
%-Strip off blank line from str2mat concatenations
if size(Gc,2)
Gcnames(1,:)=[];
end
end
%-Since no FxC interactions these are the same
Gnames = Gcnames;
%-Work out exchangability blocks
%-----------------------------------------------------------------------
%-Valid sizes of X-blk (Xblk) are integer divisors of nScan that are >1
tmp = nScan./[nScan:-1:1];
Xblk = tmp( floor(tmp)==ceil(tmp) & tmp>1 );
tmp = int2str(Xblk(1));
for i=2:length(Xblk), tmp=str2mat(tmp,int2str(Xblk(i))); end
if length(Xblk)>1
Xblk = job.xblock;%spm_input('Size of exchangability block','+0','b',tmp,Xblk);
end
nXblk = (nScan/Xblk);
iXblk = meshgrid(1:nXblk,1:Xblk); iXblk = iXblk(:)';
%-Work out how many perms, and ask about approximate tests
%-----------------------------------------------------------------------
%-NB: n! == gamma(n+1)
nPiCond_mx = round(exp(nXblk*gammaln(Xblk+1)));
% bAproxTst = spm_input(sprintf('%d Perms. Use approx. test?',nPiCond),...
% '+1','y/n')=='y';
nPiCond = job.nPerm;
if job.nPerm >= nPiCond_mx
bAproxTst=0;
if job.nPerm > nPiCond_mx
fprintf('NOTE: %d permutations requested, only %d possible.\n',job.nPerm, nPiCond_mx)
nPiCond = nPiCond_mx;
end
else
bAproxTst=1;
end
%-Compute permutations of conditions
%=======================================================================
if bAproxTst
%-Approximate test :
% Build up random subset of all (within Xblk) permutations
%===============================================================
PiCond = zeros(nPiCond,nScan);
PiCond(1,:) = 1+rem([0:Xblk*nXblk-1],Xblk);
for i = 2:nPiCond
%-Generate a new random permutation - see randperm
[null,p] = sort(rand(Xblk,nXblk)); p = p(:)';
%-Check it's not already in PiCond
while any(all((meshgrid(p,1:i-1)==PiCond(1:i-1,:))'))
[null,p] = sort(rand(Xblk,nXblk)); p = p(:)';
end
PiCond(i,:) = p;
end
clear p
else
%-Full permutation test :
% Build up exhaustive matrix of permutations
%===============================================================
%-Compute permutations for a single exchangability block
%---------------------------------------------------------------
%-Initialise XblkPiCond & remaining numbers
XblkPiCond = [];
lef = [1:Xblk]';
%-Loop through numbers left to add to permutations, accumulating PiCond
for i = Xblk:-1:1
%-Expand XblkPiCond & lef
tmp = round(exp(gammaln(Xblk+1)-gammaln(i+1)));
Exp = meshgrid(1:tmp,1:i); Exp = Exp(:)';
if ~isempty(XblkPiCond), XblkPiCond = XblkPiCond(:,Exp); end
lef = lef(:,Exp);
%-Work out sampling for lef
tmp1 = round(exp(gammaln(Xblk+1)-gammaln(i+1)));
tmp2 = round(exp(gammaln(Xblk+1)-gammaln(i)));
sam = 1+rem(0:i*tmp1-1,i) + ([1:tmp2]-1)*i;
%-Add samplings from lef to XblkPiCond
XblkPiCond = [XblkPiCond; lef(sam)];
%-Delete sampled items from lef & condition size
lef(sam) = [];
tmp = round(exp(gammaln(Xblk+1)-gammaln((i-1)+1)));
lef = reshape(lef,(i-1),tmp);
%NB:gamma(Xblk+1)/gamma((i-1)+1) == size(XblkPiCond,2);
end
clear lef Exp sam i
%-Reorientate so permutations are in rows
XblkPiCond = XblkPiCond';
%-Now build up the complete set of permutations: PiConds
%---------------------------------------------------------------
nXblkPiCond=size(XblkPiCond,1);
PiCond=XblkPiCond;
for i=2:nXblk
tmp=size(PiCond,1);
[iSup,iInf] = meshgrid(1:nXblkPiCond,1:tmp);
iSup = iSup(:)';
iInf = iInf(:)';
PiCond=[XblkPiCond(iSup,:),PiCond(iInf,:)];
end
end
%-Check, condition and randomise PiCond
%-----------------------------------------------------------------------
%-Check PiConds sum within Xblks to sum of first nXblk natural numbers
if ~all(all(PiCond*spm_DesMtx(iXblk)== (Xblk+1)*Xblk/2 ))
error('SnPM:InvalidPiCond', 'Invalid PiCond computed!'), end
%-Convert to full permutations from permutations within blocks
nPiCond = size(PiCond,1);
PiCond = PiCond + meshgrid((iXblk-1)*Xblk,1:nPiCond);
%-Randomise order of PiConds (except first) to allow interim analysis
PiCond=[PiCond(1,:);PiCond(randperm(nPiCond-1)+1,:)];
%-Check first permutation is null permutation
if ~all(PiCond(1,:)==[1:nScan])
error('SnPM:InvalidPiCond', 'PiCond(1,:)~=[1:nScan]'); end
%-Form non-null design matrix partitions (Globals handled later)
%=======================================================================
CONT = [1]; % Contrast of condition effects
%-Covariate partition & Form for HC computation at permutation perm
if bCCovInt
C = CovInt - mean(CovInt);
Cnames = 'CovInt(centred)';
Cc = CovInt;
Ccnames = 'CovInt';
sHCform = ['spm_DesMtx(CovInt(PiCond(perm,:))-mean(CovInt),',...
'''C'',''CovInt(centred)'')'];
else
C = CovInt;
Cnames = 'CovInt';
sHCform = 'spm_DesMtx(CovInt(PiCond(perm,:)),''C'',''CovInt'')';
end
%-Include constant term in block partition
B=ones(nScan,1); Bnames='Const';
%-Design description
%-----------------------------------------------------------------------
sDesign = sprintf('SingleSub: Simple Regression (correlation); single covariate of interest: %d scans',nScan);
sPiCond = sprintf('Permutations of covariate within exchangability blocks of size %d, %d permutations',Xblk,size(PiCond,1));