ensemble_response_times

 Analyzes response time values stored in the response_text field of a
 response table.

 [out_st] = ensemble_response_times(data_st,params);

 The questions for which response times are to be analyzed should be
 specified in a filtering directive stored in the params structure, e.g.
 filt.include.all.question_id = [question_ids];

 params.rt.rt_fld = name of field containing response time data (in msec).
 Default is response_text

 params.rt.multiRespCheckFld = name of field to use for checking
 uniqueness of response to each stimulus, i.e. one response per stimulus
 per subject. Default='session_id';

 params.rt.transform = {'none','log10'}; - specifies whether a transform
 should be applied to the RT data. Default='none'.

 params.rt.proc = {}; - list of functions to process. The default behavior is to
 treat the name as a function handle. Default processing includes
 calculation of the mean and standard deviation.

 params.rt.by_item = {'stimulus_id','subject_id'};

 NOTE: A cleaner implementation of this would require that data are
 already pre-filtered. What is required for that to happen though is a
 modification of ensemble_filter() selectively remove rows based on
 combinations of variable values.

0001 function [out_st] = ensemble_response_times(data_st,params)
0002 % Analyzes response time values stored in the response_text field of a
0003 % response table.
0004 %
0005 % [out_st] = ensemble_response_times(data_st,params);
0006 %
0007 % The questions for which response times are to be analyzed should be
0008 % specified in a filtering directive stored in the params structure, e.g.
0009 % filt.include.all.question_id = [question_ids];
0010 %
0011 % params.rt.rt_fld = name of field containing response time data (in msec).
0012 % Default is response_text
0013 %
0014 % params.rt.multiRespCheckFld = name of field to use for checking
0015 % uniqueness of response to each stimulus, i.e. one response per stimulus
0016 % per subject. Default='session_id';
0017 %
0018 % params.rt.transform = {'none','log10'}; - specifies whether a transform
0019 % should be applied to the RT data. Default='none'.
0020 %
0021 % params.rt.proc = {}; - list of functions to process. The default behavior is to
0022 % treat the name as a function handle. Default processing includes
0023 % calculation of the mean and standard deviation.
0024 %
0025 % params.rt.by_item = {'stimulus_id','subject_id'};
0026 %
0027 % NOTE: A cleaner implementation of this would require that data are
0028 % already pre-filtered. What is required for that to happen though is a
0029 % modification of ensemble_filter() selectively remove rows based on
0030 % combinations of variable values.
0031 
0032 
0033 % 18Nov2012 Petr Janata
0034 % 05Feb2013 PJ - added flexible specification of field containing response
0035 %                times
0036 
0037 out_st = ensemble_init_data_struct;
0038 out_st.type = 'response_time_data';
0039 out_st.vars = [{'data_transform','rank'} params.rt.calc];
0040 
0041 % Get the idx of the response_data
0042 ridx = ensemble_find_analysis_struct(data_st,struct('type','response_data'));
0043 if isempty(ridx)
0044   ridx = ensemble_find_analysis_struct(data_st,struct('type','resp_x_stim'));
0045 end
0046 if isempty(ridx)
0047   error('Could not find an appropriate data structure with response data');
0048 end
0049 rst = data_st{ridx};
0050 rcols = set_var_col_const(rst.vars);
0051 
0052 % Apply an standard filtering parameters
0053 if isfield(params,'filt')
0054   fprintf('Applying filtering criteria\n')
0055   rst = ensemble_filter(rst, params.filt);
0056 end
0057 
0058 % See if there is a datastruct containing filtering parameters and apply
0059 % them
0060 if isfield(params.rt,'exclude')
0061   rm_vars = {'name','type','report','meta'};
0062   idx = ensemble_find_analysis_struct(data_st,struct('type','filtering_data'));
0063   
0064   fparams = rmfield(ensemble_tree2datastruct(data_st{idx}), rm_vars);
0065   nfilt = size(fparams.data{1},1);
0066   
0067   % Have to filter row by row
0068   for ifilt = 1:nfilt
0069     tmp.vars = fparams.vars;
0070     for ivar = 1:length(fparams.vars)
0071       tmp.data{ivar} = fparams.data{ivar}(ifilt);
0072     end
0073     f.exclude.all = rmfield(ensemble_datastruct2tree(tmp), rm_vars);
0074     rst = ensemble_filter(rst,f);
0075   end
0076 end
0077 
0078 % Determine whether we are looping by a particular item (stimulus_id,
0079 % subject_id)
0080 if isfield(params,'item_var')
0081   item_str = params.item_var;
0082 else
0083   item_str = 'stimulus_id';
0084 end
0085 
0086 switch item_str
0087   case 'session_id'
0088     item_type = 'session';
0089   case 'subject_id'
0090     item_type = 'subject';
0091   case 'stimulus_id'
0092     item_type = 'stimulus';
0093   case 'none'
0094     item_type = 'none';
0095   otherwise
0096     fprintf('%s: Unknown item variable: %s\n', mfilename, item_str);
0097 end
0098 
0099 % Set the output structure name
0100 if isfield(params, 'outputNameType')
0101   outputNameType = params.outputNameType;
0102 else
0103   outputNameType = 'itemType';
0104 end
0105 
0106 switch outputNameType
0107   case 'itemType'
0108     out_st.name = sprintf('by_%s', item_str);
0109   case 'analysisName'
0110     out_st.name = '';  % allow ensemble_jobman to populate it
0111   otherwise
0112     out_st.name = '';
0113 end
0114 
0115 % Add item types to the list of output variables
0116 out_st.vars = [item_str out_st.vars];
0117 ocols = set_var_col_const(out_st.vars);
0118 
0119 % If we are looping over stimulus_id and we are normalizing the response
0120 % time data by some value, make sure we have a structure of those values
0121 if strcmp(item_type,'stimulus') && isfield(params.rt,'transform') && any(ismember(params.rt.transform,'norm'))
0122   % Get the stimulus info
0123   sidx = ensemble_find_analysis_struct(data_st,struct('type','stimulus_metadata'));
0124   if isempty(sidx)
0125     error('stimulus_metadata structure required for analysis, but none found\n')
0126   end
0127   stim_st = data_st{sidx};
0128   stimcols = set_var_col_const(stim_st.vars);
0129 end
0130 
0131 % Precalculate the masks we are using for this looping variable
0132 if strcmp(item_type,'none')
0133   item_mask_mtx = ones(length(rst.data{1}),1);
0134   itemids = 1;
0135 else
0136   % Precalculate the masks for the item variable we are using
0137   [item_mask_mtx, itemids] = make_mask_mtx(rst.data{rcols.(item_str)});
0138 end
0139 nitems = length(itemids);
0140 
0141 % Determine which field the response times are in
0142 if isfield(params.rt, 'rt_fld')
0143   rt_fld = params.rt.rt_fld;
0144 else
0145   rt_fld = 'response_text';
0146 end
0147 
0148 % Convert the response text data to doubles
0149 if iscell(rst.data{rcols.(rt_fld)})
0150   allRTs = cellfun(@str2num,rst.data{rcols.(rt_fld)});
0151 else
0152   allRTs = rst.data{rcols.(rt_fld)};
0153 end
0154 
0155 % Convert to seconds
0156 allRTs = allRTs/1000;
0157 
0158 if isfield(params.rt, 'multiRespCheckFld')
0159   multiRespCheckFld = params.rt.multiRespCheckFld;
0160 else
0161   multiRespCheckFld = 'session_id';
0162 end
0163 
0164 % Now loop over items
0165 for iitem = 1:nitems
0166   item_mask = item_mask_mtx(:,iitem);
0167   
0168   % Check to see if we have multiple responses per subject per item
0169   uniqueSess = unique(rst.data{rcols.(multiRespCheckFld)}(item_mask));
0170   if iscell(uniqueSess)
0171     [~,cnt] = cellhist(uniqueSess);
0172   else
0173     [cnt] = hist(rst.data{rcols.(multiRespCheckFld)}(item_mask),uniqueSess);
0174   end
0175   
0176   if any(cnt > 1)
0177     multmask = ismember(rst.data{rcols.(multiRespCheckFld)}, uniqueSess(cnt > 1));
0178     uniqueMult = unique(rst.data{rcols.(multiRespCheckFld)}(multmask));
0179     nmult = length(uniqueMult);
0180     for imult = 1:nmult
0181       currMult = uniqueMult(imult);
0182       currmask = ismember(rst.data{rcols.(multiRespCheckFld)}, currMult) & item_mask;
0183 
0184       fprintf('Found %d responses for stimulus %d, session %d\n', ...
0185         sum(currmask), itemids(iitem), currMult);
0186 
0187       % Take the first of the multiple responses
0188       [~,minidx] = min(rst.data{rcols.date_time}(currmask));
0189       curridxs = find(currmask);
0190       goodmask = false(size(currmask));
0191       goodmask(curridxs(minidx)) = true;
0192       item_mask = xor(item_mask,xor(currmask,goodmask));
0193     end
0194   end
0195   
0196   % Extract the response times
0197   resptimes = allRTs(item_mask);
0198    
0199   % Determine whether we are going to transform the data
0200   if isfield(params.rt,'transform')
0201     xfms = params.rt.transform;
0202     nxfm = length(xfms);
0203   else
0204     nxfm = 1;
0205     xfms = {'none'};
0206   end
0207   
0208   % Loop over possible transformations of the data
0209   for ixfm = 1:nxfm
0210     xfmType = xfms{ixfm};
0211     
0212     row = (iitem-1)*nxfm + ixfm; % Determine which row we are putting data into
0213     
0214     % Output information about the item ID
0215     if iscell(itemids)
0216       out_st.data{ocols.(item_str)}{row,1} = itemids(iitem);
0217     else
0218       out_st.data{ocols.(item_str)}(row,1) = itemids(iitem);
0219     end
0220     
0221     % Output information about the transform type
0222     out_st.data{ocols.data_transform}{row,1} = xfmType;
0223     
0224     % Transform the data
0225     switch xfmType
0226       case 'none'
0227         data = resptimes;
0228       case 'norm' % normalize by item-level data
0229         % Normalize by a particular field
0230         if strcmp(params.rt.norm,'duration') && isnan(stim_st.data{stimcols.duration}(iitem))
0231           srcmask = stim_st.data{stimcols.stimulus_id} == itemids(iitem);
0232           stimname = fullfile(params.paths.stimulus_root, stim_st.data{stimcols.location}{srcmask});
0233           system_str = sprintf('mp3info -p %%S %s', stimname);
0234           [~,r] = system(system_str);
0235           stim_st.data{stimcols.duration}(srcmask) = min(params.maxListenTime,str2double(r));
0236         end
0237 
0238         normVal = stim_st.data{stimcols.(params.rt.norm)}(stim_st.data{stimcols.stimulus_id} == itemids(iitem));
0239         data = resptimes/normVal;
0240       otherwise
0241         fh = str2func(xfmType);
0242         data = fh(resptimes);
0243     end
0244     
0245     % Calculate the desired metrics on these response times
0246     ncalc = length(params.rt.calc);
0247     for icalc = 1:ncalc
0248       currCalc = params.rt.calc{icalc};
0249       
0250       switch currCalc
0251         
0252         otherwise
0253           fh = str2func(currCalc);
0254           out_st.data{ocols.(currCalc)}(row,1) = fh(data);
0255           
0256       end % switch currCalc
0257     end % for icalc
0258   end % for ixfm
0259 end % for iitem
0260 
0261 % Perform additional analyses across items by transform type
0262 rankVars = {'numel','std',item_str};
0263 for ixfm = 1:nxfm
0264   xfmType = xfms{ixfm};
0265   xfm_mask = ismember(out_st.data{ocols.data_transform}, xfmType);
0266   
0267   % Get the rank order of the mean data
0268   [ranked(ixfm).mean, srcidxs] = sort(out_st.data{ocols.mean}(xfm_mask),'descend');
0269   
0270   % Write the rank info to the by-item data matrix
0271   [sorted_srcidx, ranks] = sort(srcidxs);
0272   out_st.data{ocols.rank}(xfm_mask,1) = ranks; 
0273   
0274   for ivar = 1:length(rankVars)
0275     cv = rankVars{ivar};
0276     
0277     % Rank the current variable
0278     tmp = out_st.data{ocols.(cv)}(xfm_mask);
0279     ranked(ixfm).(cv) = tmp(srcidxs);
0280   end
0281 
0282   % Calculate the standard error of the mean
0283   ranked(ixfm).sem = ranked(ixfm).std./sqrt(ranked(ixfm).numel-1);
0284   
0285   % Generate a rank ordered bar graph
0286   subplot(nxfm,1,ixfm)
0287   h = bar(ranked(ixfm).mean,'facecolor',ones(1,3)*0.8);
0288   add_errorbars(h,ranked(ixfm).sem,'k')
0289   
0290   % Add xtick labels
0291   set(gca,'xtick',1:nitems,'xticklabel','')
0292   
0293   for iitem = 1:nitems
0294     if iscell(ranked(ixfm).(item_str))
0295       itemStr = sprintf('%s',ranked(ixfm).(item_str){iitem});
0296     else
0297       itemStr = sprintf('%d',ranked(ixfm).(item_str)(iitem));
0298     end
0299     text(iitem,max(get(gca,'ylim'))*0.05,itemStr,'rotation',90)
0300   end
0301   
0302   xlabel(strrep(item_str,'_','\_'), 'fontsize', 14, 'fontweight', 'bold')
0303   
0304   switch lower(xfmType)
0305     case 'none'
0306       lblstr = 'Response time (s)';
0307     otherwise
0308       lblstr = sprintf('%s(s)', xfmType);
0309   end
0310   ylabel(lblstr, 'fontsize', 14, 'fontweight','bold')
0311   
0312   % Add sample size info
0313   if ~any(diff(ranked(ixfm).numel))
0314     nobsstr = sprintf('N = %d',ranked(ixfm).numel(1));
0315   else
0316     nobsstr = sprintf('N = %d - %d', min(ranked(ixfm).numel), max(ranked(ixfm).numel));
0317   end
0318   text(0.95,0.95,nobsstr,'units','normalized',...
0319     'horizontalalign','right', ...
0320     'fontsize', 16, ...
0321     'fontweight','bold')
0322 end
0323 
0324 end % function
0325