mulitple regression in 2nd level batch

Author: torwager

Second_level_analysis_template_scripts/b_copy_to_local_scripts_dir_and_modify/a2_set_default_options.m

@@ -14,6 +14,9 @@
 % --------------------------------------------------------------------
 dorobust = true;            % robust statistics [true, false] -- default true
 myscaling = 'raw';          % 'raw' or 'scaled'
+design_matrix_type = 'custom';   % 'group' or 'custom'
+                            % Group: use DAT.BETWEENPERSON.group or DAT.BETWEENPERSON.contrasts{c}.group;
+                            % Custom: use all columns of table object DAT.BETWEENPERSON.contrasts{c};
 
 % prep_3b_run_SVMs_on_contrasts_and_save options 
 % --------------------------------------------------------------------

Second_level_analysis_template_scripts/core_scripts_to_run_without_modifying/d9_nps_correlations_with_global_signal.m

@@ -9,9 +9,28 @@
 
     clear r p
 
+    % corr performance changed from Matlab 2017b to 2018a. Could be bug in
+    % 2018a.  
+    use_alt_method = false;
+    myversion = ver('Matlab'); 
+    if datenum(myversion.Date) > datenum('01-Jan-2018')
+        use_alt_method = true;
+    end
+    
+    
     for i = 1:k
 
-        [r(i, :), p(i, :)] = corr(DAT.npscontrasts{i}, DAT.gray_white_csf_contrasts{i});
+        if use_alt_method
+            
+            [tmpr, tmpp] = corr([DAT.npscontrasts{i} DAT.gray_white_csf_contrasts{i}]);
+            r(i, :) = tmpr(1, 2:end);
+            p(i, :) = tmpp(1, 2:end);
+            
+        else
+            
+            [r(i, :), p(i, :)] = corr(DAT.npscontrasts{i}, DAT.gray_white_csf_contrasts{i});
+            
+        end
 
     end
 

Second_level_analysis_template_scripts/core_scripts_to_run_without_modifying/prep_3a_run_second_level_regression_and_save.m

@@ -1,12 +1,29 @@
-% THIS SCRIPT RUNS BETWEEN-PERSON CONTRASTS using second-level regression
+% THIS SCRIPT RUNS BETWEEN-PERSON (2nd-level) Regression analyses
+% for each within-person CONTRAST registered in the analysis
+% 
+% - To specify analysis options, run a2_set_default_options
+% - prep_3a_run_second_level_regression_and_save runs regressions and saves
+% results in a standard location and format
+% - To get results reports, see c2a_second_level_regression
 %
-% FOR NOW: Uses "group" variable.
-% FUTURE: Extend this to arbitrary covariates (> 1 cov), which have to be
-% entered in DAT structure in some standardize place.
+% Analysis options include:
+% - dorobust : robust regression or OLS (true/false)
+% - myscaling: 'raw' or 'scaled' (image scaling done in prep_2_... data load)
+% - design_matrix_type: 'group' or 'custom'
+%                       Group: use DAT.BETWEENPERSON.group or DAT.BETWEENPERSON.contrasts{c}.group;
+%                       Custom: use all columns of table object DAT.BETWEENPERSON.contrasts{c};
 %
-% Assuming that groups are concatenated into contrast image lists.
-% Requires DAT.BETWEENPERSON.group field specifying group membership for
+% 'group' option 
+% Assuming that groups are concatenated in contrast image lists, and
+% regressor values of 1 or -1 will specify the group identity for each
+% image. Requires DAT.BETWEENPERSON.group field specifying group membership for
 % each image.
+%
+% 'custom' option: 
+% Can enter a multi-column design matrix for each contrast
+% Design matrix can be different for each contrast
+%
+% To set up group and custom variables, see prep_1b_prep_behavioral_data
 
 % --------------------------------------------------------------------
 
@@ -21,10 +38,10 @@
 % (e.g., not specified in an older/incomplete copy of a2_set_default_options)
 
 % Now set in a2_set_default_options
-options_needed = {'dorobust', 'myscaling'};  % Options we are looking for. Set in a2_set_default_options
+options_needed = {'dorobust', 'myscaling', 'design_matrix_type'};  % Options we are looking for. Set in a2_set_default_options
 options_exist = cellfun(@exist, options_needed); 
 
-option_default_values = {false 'raw'};          % defaults if we cannot find info in a2_set_default_options at all 
+option_default_values = {false, 'raw', 'group'};          % defaults if we cannot find info in a2_set_default_options at all 
 
 plugin_get_options_for_analysis_script
 
@@ -65,15 +82,34 @@
 
 for c = 1:kc
 
+    % Get design matrix for this contrast
+    % --------------------------------------------------------------------
+    
     mygroupnamefield = 'contrasts';  % 'conditions' or 'contrasts'
-    [group, groupnames, groupcolors] = plugin_get_group_names_colors(DAT, mygroupnamefield, c);
-    outcome_value = group;
 
-    if isempty(group)
-        fprintf('Group not defined for contrast %s. Skipping.\n', DAT.contrastnames{c}); 
-        continue
+    switch design_matrix_type
+        case 'custom'
+            
+            % Define design matrix X "design_matrix"
+            % Use custom matrix for each condition/contrast
+            table_obj = DAT.BETWEENPERSON.(mygroupnamefield){c};
+            groupnames = table_obj.Properties.VariableNames;
+            X = table2array(table_obj);
+            
+        case 'group'
+            
+            % Use 'groups' single regressor
+            [group, groupnames, groupcolors] = plugin_get_group_names_colors(DAT, mygroupnamefield, c);
+            X = group;
+            
+            if isempty(group)
+                fprintf('Group not defined for contrast %s. Skipping.\n', DAT.contrastnames{c});
+                continue
+            end
+            
+        otherwise error('Incorrect option specified for design_matrix_type');
     end
-
+    
     printstr(DAT.contrastnames{c});
     printstr(dashes)
 
@@ -100,26 +136,58 @@
 
     % a. Format and attach outcomes: 1, -1 group variable FOR NOW
     % --------------------------------------------------------------------
- 
-    % Confirm outcome_value is 1, -1, or mean-centered
-    % *****add this*****
 
-    % Check multilinearity, variance in all regressors; error for bad
-    % design matrix
-    % *****add this*****
+    % Confirm design_matrix is 1, -1, or mean-centered
+    meancentered = ~(abs(mean(X)) > 1000 * eps);
+    effectscoded = all(X == 1 | X == -1 | X == 0, 1);
+    isconstant = all(X == mean(X, 1), 1);
+    vifs = getvif(X);
 
-    cat_obj.X = outcome_value;
+    if any(isconstant)
+        disp('An intercept appears to be added manually. Do not include an intercept - it will be added automatically.');
+        disp('Skipping this contrast.')
+        continue
+    end
+
+    % Report
+    design_table = table;
+    design_table.Mean = mean(X)';
+    design_table.Var = var(X)';
+    design_table.EffectsCode = effectscoded';
+    design_table.VIF = vifs';
+    design_table.Properties.RowNames = groupnames';
+    disp(design_table)
+    disp(' ');
 
-    % Skip if necessary
-    % --------------------------------------------------------------------
+    if any(~meancentered & ~effectscoded)
+        disp('Warning: some columns are not mean-centered or effects coded. \nIntercept may not be interpretable.\n');
+        fprintf('Columns: ')
+        fprintf('%d ', find(~meancentered & ~effectscoded));
+        fprintf('\n');
+    else
+        disp('Checked OK: All columns mean-centered or are effects-coded [1 -1 0]');
+    end
 
-    if all(cat_obj.X > 0) || all(cat_obj.X < 0)
-        % Only positive or negative weights - nothing to compare
-        
-        printhdr(' Only positive or negative regressor values - bad design');
-        
-        continue
+    
+    if any(vifs) > 2
+        disp('Some regressors have high variance inflation factors: Parameters might be poorly estimated or uninterpretable.');
+    else
+        disp('Checked OK: VIFs for all columns are < 2');
     end
+
+    
+    cat_obj.X = X;
+    
+    % Skip if necessary
+    % --------------------------------------------------------------------
+%     
+%     if all(cat_obj.X > 0) || all(cat_obj.X < 0)
+%         % Only positive or negative weights - nothing to compare
+%         
+%         printhdr(' Only positive or negative regressor values - bad design');
+%         
+%         continue
+%     end
 
     % Run voxel-wise regression model
     % --------------------------------------------------------------------
@@ -134,7 +202,9 @@
     % out.t has t maps for all regressors, intercept is last
     regression_stats = regress(cat_obj, .05, 'unc', robuststring, 'analysis_name', DAT.contrastnames{c}, 'variable_names', {'Group'});
 
+
     % Make sure variable types are right data formats
+    regression_stats.design_table = design_table;
     regression_stats.t = enforce_variable_types(regression_stats.t);
     regression_stats.b = enforce_variable_types(regression_stats.b);
     regression_stats.df = enforce_variable_types(regression_stats.df);
@@ -146,7 +216,12 @@
 
     % add names for analyses and variables: 
     regression_stats.analysis_name = DAT.contrastnames{c};
-    regression_stats.variable_names = {'Group' 'Intercept'}
+    regression_stats.variable_names = [groupnames {'Intercept'}];
+    
+    % prints output automatically - name axis
+    for kk = 1:length(regression_stats.variable_names)
+        spm_orthviews_name_axis(regression_stats.variable_names{kk}, kk);
+    end
 
     % Save stats objects for results later
     % --------------------------------------------------------------------
@@ -164,6 +239,7 @@
 
 save(savefilenamedata, 'regression_stats_results', '-v7.3');
 printhdr('Saved regression_stats_results for contrasts');
+fprintf('Filename: %s\n', savefilenamedata);