simulate_model_norm

 prime_base_model_norm
 
 This script sets up an event-related analysis.  Although the model for each
 subject is set up from the same prototype, models will differ across
 subjects. One MAJOR DIFFERENCE in models across subjects may be the number
 of conditions modeled.  Since all targets are modeled as a function of
 response accuracy, if subjects made no errors for a particular type of target
 in a particular session, the column modeling that condition will be removed
 from the design matrix.  This means not only that the d.f. will differ across
 subjects but also that the contrasts have to be set up on a subject by
 subject basis.

0001 % prime_base_model_norm
0002 %
0003 % This script sets up an event-related analysis.  Although the model for each
0004 % subject is set up from the same prototype, models will differ across
0005 % subjects. One MAJOR DIFFERENCE in models across subjects may be the number
0006 % of conditions modeled.  Since all targets are modeled as a function of
0007 % response accuracy, if subjects made no errors for a particular type of target
0008 % in a particular session, the column modeling that condition will be removed
0009 % from the design matrix.  This means not only that the d.f. will differ across
0010 % subjects but also that the contrasts have to be set up on a subject by
0011 % subject basis.
0012 %
0013 
0014 
0015 global dataroot spm_path subject_ids use_model_proto sinfo
0016 global SCAN_OFFSET
0017 global REGRESS_MOTION_CORRECTION_PARAMS ADD_LINEAR_TREND_REGRESS ADD_RT_REGRESS
0018 
0019 nsub = length(subject_ids);
0020 
0021 subj_nsess = length(sinfo(1).cond_order);
0022 
0023 disp('Getting subject file info for models')
0024 for isub = 1:nsub
0025   for isess = 1:subj_nsess
0026     sessions(isub).images{isess} = spm_get('Files', ... 
0027     sprintf('%s/%s/epi/run%d/', dataroot, char(subject_ids(isub)), isess), ...
0028     sprintf('sn%s*.img', char(subject_ids(isub))) ...
0029     );   
0030   end  % for isess=
0031 end  % for isub=
0032 
0033 prime_long_model_protos;
0034 
0035 model_init = model_prototype(use_model_proto);
0036 conditions_init = conditions_prototype(use_model_proto);
0037  
0038 if REGRESS_MOTION_CORRECTION_PARAMS | ADD_LINEAR_TREND_REGRESS
0039   regressors_init = regressors_prototype(use_model_proto);
0040 end
0041 
0042 ncond_idxs = 0;
0043 nregress_idxs = 0;
0044   
0045 clear regressors
0046 
0047 disp('Getting subject specific model parameters')
0048 for isub = 1:nsub
0049   model(isub) = model_init;
0050 
0051   model(isub).files = sessions(isub).images;
0052   
0053   % Set the SOT values for this subject
0054   ttinfo_fname = fullfile(dataroot, char(subject_ids(isub)), sprintf('%s_ttinfo.mat', char(subject_ids(isub))));
0055   sot = prime_SOTs(ttinfo_fname, use_model_proto);
0056   
0057   for isess = 1:model(isub).nsess
0058     ncond_idxs = ncond_idxs + 1;
0059     sub_conds(isess) = ncond_idxs;        % Keep track of which condition numbers
0060                                         % correspond to the model for this subject
0061 
0062     conditions(ncond_idxs) = conditions_init;
0063     conditions(ncond_idxs).onsets = sot{isess};
0064     
0065     %
0066     % Check to make sure that there are events in every condition, and remove
0067     % conditions for which there are no events
0068     %
0069     
0070     % Compile list of empty conditions
0071     ncond = length(sot{isess});
0072     empty_cond = [];
0073     for icond = 1:ncond
0074       if isempty(sot{isess}{icond})
0075     empty_cond(end+1) = icond;
0076       end
0077     end
0078     
0079     %
0080     % Delete entries from cond_names, cond_types, onsets, bf_ev, and bf_ep
0081     %
0082     
0083     if empty_cond
0084       conditions(ncond_idxs).names(empty_cond) = [];
0085       conditions(ncond_idxs).onsets(empty_cond) = [];
0086       conditions(ncond_idxs).types(empty_cond) = [];
0087       conditions(ncond_idxs).bf_ev(empty_cond) = [];
0088       conditions(ncond_idxs).bf_ep(empty_cond) = [];
0089       conditions(ncond_idxs).durations{1}(empty_cond) = [];
0090     end % if empty_cond
0091     
0092     sub_nconds(isess) = length(conditions(ncond_idxs).names);
0093 
0094     if REGRESS_MOTION_CORRECTION_PARAMS | ADD_LINEAR_TREND_REGRESS
0095       nregress_idxs = nregress_idxs + 1;
0096       regressors(nregress_idxs) = regressors_init;
0097       nregress_sess = 0;
0098     end
0099     
0100       
0101     if REGRESS_MOTION_CORRECTION_PARAMS
0102       % Load regressors
0103       %fname = spm_get('Files',sprintf('%s/%s/epi/run%d/', dataroot, char(subject_ids(isub)), isess),sprintf('realignment_params_*i%04d.txt', SCAN_OFFSET+1));
0104       fname = spm_get('Files',sprintf('%s/%s/epi/run%d/', dataroot, char(subject_ids(isub)), isess),sprintf('realignment_params_%s*i%04d.txt', char(subject_ids(isub)),SCAN_OFFSET+1));
0105       vals = load(fname);
0106       ncol = size(vals,2);
0107       regressors(nregress_idxs).values(1:sinfo(isub).nvol(isess)-SCAN_OFFSET,nregress_sess+1:nregress_sess+ncol) = load(fname);
0108       nregress_sess = nregress_sess + ncol;
0109     end
0110     
0111     if ADD_LINEAR_TREND_REGRESS
0112       vals = 0:sinfo(isub).nvol(isess)-SCAN_OFFSET-1;
0113       vals = vals'/max(vals);
0114       regressors(nregress_idxs).values(1:sinfo(isub).nvol(isess)-SCAN_OFFSET,nregress_sess+1) = vals;
0115       nregress_sess = nregress_sess + 1;
0116     end
0117     
0118     if ADD_RT_REGRESS
0119       % load RTs for regressor
0120       a_info=load(ttinfo_fname); 
0121       % how many RTs per run?
0122       nb = length(a_info.rt.raw)/model(isub).nsess;
0123       a_start = (isess-1)*nb+1; a_stop = isess*nb;
0124       % for rcr
0125       a_idx_rcr = find(a_info.good_idx.rcr >= a_start & a_info.good_idx.rcr <= a_stop);
0126       a_rcr = a_info.good_idx.rcr(a_idx_rcr);
0127       a_help = a_info.rt.raw(a_rcr);
0128       a_when = round(sot{isess}{4});
0129       vals(a_when) = a_help; clear a_help a_when
0130       % for ucr
0131       a_idx_ucr = find(a_info.good_idx.ucr >= a_start & a_info.good_idx.ucr <= a_stop);
0132       a_ucr = a_info.good_idx.ucr(a_idx_ucr);
0133       a_help = a_info.rt.raw(a_ucr);
0134       a_when = round(sot{isess}{6});
0135       vals(a_when) = a_help;clear a_help a_when
0136       % for rdr
0137       a_idx_rdr = find(a_info.good_idx.rdr >= a_start & a_info.good_idx.rdr <= a_stop);
0138       a_rdr = a_info.good_idx.rdr(a_idx_rdr);
0139       a_help = a_info.rt.raw(a_rdr);
0140       a_when = round(sot{isess}{5});
0141       vals(a_when) = a_help;clear a_help a_when
0142       % for udr
0143       a_idx_udr = find(a_info.good_idx.udr >= a_start & a_info.good_idx.udr <= a_stop);
0144       a_udr = a_info.good_idx.udr(a_idx_udr);
0145       a_help = a_info.rt.raw(a_udr);
0146       a_when = round(sot{isess}{7});
0147       vals(a_when) = a_help;clear a_help a_when
0148 
0149       % eliminate the NaN of outliers
0150       a_nan = isnan(vals); a_wo = find(a_nan > 0); 
0151       for ai=1:length(a_wo)
0152     vals(a_wo(ai,1)) = 0;  % eliminate the NaN of outliers
0153       end 
0154       
0155       regressors(nregress_idxs).values(1:sinfo(isub).nvol(isess)-SCAN_OFFSET,nregress_sess+1) = vals;
0156       nregress_sess = nregress_sess + 1;
0157     end  
0158 
0159     sub_regress(isess) = nregress_idxs;
0160     sub_nregress(isess) = length(regressors(nregress_idxs).names);
0161   end % for isess =
0162   
0163   % Set the condition info for this subject
0164   model(isub).conditions = sub_conds;
0165   model(isub).conditions_nb = sub_nconds;
0166   
0167   % Set the regressor info for this subject
0168   model(isub).regressors = sub_regress;  
0169   model(isub).regressors_nb = sub_nregress;
0170   
0171   % Set the correct number of scans for this subject
0172   model(isub).nscans = sinfo(isub).nvol - SCAN_OFFSET;
0173   
0174 end  %for isub
0175