% LMM code for Hickey, Grignolio, Acunzo, & Munasinghe (2024) % requires MATLAB statistics toolbox tbl = readtable('table.csv'); tbl.condition = cell2mat(tbl.condition); load_models_from_file = true; if load_models_from_file load('models.mat'); else %%% Combined analysis of target and distractor Nt % The random effects have a nested structure. This is equivalent to % analysis with independent labels for each of the 160 images. Object % category is not treated as an independent random effect as it has too few % levels (consensus is random effects need >6 levels). Category % can be included but this does not improve models, factors get removed % during parameter selection. % Note that modelling is (implicitly) using DummyVarCoding 'reference' % and this is important in interpreting results from the omnibus. % The reference level is the 'targ' condition. The test of 'RT' is % conducted in this condition (ignoring results from the 'dist' % condition). omnibus_model.base = fitlme(tbl, 'N2pc ~ RT * condition + (RT * condition|subjects) ', 'Verbose', 1, 'FitMethod', 'REML'); base_AIC = omnibus_model.base.ModelCriterion.AIC; omnibus_model.mdlA1 = fitlme(tbl, ['N2pc ~ RT * condition + (RT * condition|subjects) + ' ... '(RT * condition|targ_cat:targ_index) + (RT * condition|dist_cat:dist_index)'], 'Verbose', 1, 'FitMethod', 'REML'); omnibus_model.mdlA2 = fitlme(tbl, ['N2pc ~ RT * condition + (RT * condition|subjects) + ' ... '(RT * condition|targ_cat:targ_index) + (RT + condition|dist_cat:dist_index)'], 'Verbose', 1, 'FitMethod', 'REML'); omnibus_model.mdlA3 = fitlme(tbl, ['N2pc ~ RT * condition + (RT * condition|subjects) + ' ... '(RT * condition|targ_cat:targ_index) + (condition|dist_cat:dist_index)'], 'Verbose', 1, 'FitMethod', 'REML'); omnibus_model.mdlA4 = fitlme(tbl, ['N2pc ~ RT * condition + (RT * condition|subjects) + ' ... '(RT * condition|targ_cat:targ_index) + (1|dist_cat:dist_index)'], 'Verbose', 1, 'FitMethod', 'REML'); AIC_values = [omnibus_model.mdlA1.ModelCriterion.AIC omnibus_model.mdlA2.ModelCriterion.AIC ... omnibus_model.mdlA3.ModelCriterion.AIC omnibus_model.mdlA4.ModelCriterion.AIC]; mdlA_AIC = find(AIC_values == min(AIC_values)); % check if the optimal mdlA is better than the base model % mdlA_AIC - base_AIC % A3 wins omnibus_model.mdlB1 = fitlme(tbl, ['N2pc ~ RT * condition + (RT * condition|subjects) + ' ... '(RT + condition|targ_cat:targ_index) + (condition|dist_cat:dist_index)'], 'Verbose', 1, 'FitMethod', 'REML'); omnibus_model.mdlB2 = fitlme(tbl, ['N2pc ~ RT * condition + (RT * condition|subjects) + ' ... '(condition|targ_cat:targ_index) + (condition|dist_cat:dist_index)'], 'Verbose', 1, 'FitMethod', 'REML'); omnibus_model.mdlB3 = fitlme(tbl, ['N2pc ~ RT * condition + (RT * condition|subjects) + ' ... '(1|targ_cat:targ_index) + (condition|dist_cat:dist_index)'], 'Verbose', 1, 'FitMethod', 'REML'); AIC_values = [omnibus_model.mdlB1.ModelCriterion.AIC omnibus_model.mdlB2.ModelCriterion.AIC ... omnibus_model.mdlB3.ModelCriterion.AIC]; mdlB_AIC = find(AIC_values == min(AIC_values)); % check if the optimal mdlB is better than the optimal mdlA % mdlB_AIC - mdlA_AIC % B2 wins omnibus_model.mdlC1 = fitlme(tbl, ['N2pc ~ RT * condition + (RT + condition|subjects) + ' ... '(condition|targ_cat:targ_index) + (condition|dist_cat:dist_index)'], 'Verbose', 1, 'FitMethod', 'REML'); omnibus_model.mdlC2 = fitlme(tbl, ['N2pc ~ RT * condition + (condition|subjects) + ' ... '(condition|targ_cat:targ_index) + (condition|dist_cat:dist_index)'], 'Verbose', 1, 'FitMethod', 'REML'); omnibus_model.mdlC3 = fitlme(tbl, ['N2pc ~ RT * condition + (1|subjects) + ' ... '(condition|targ_cat:targ_index) + (condition|dist_cat:dist_index)'], 'Verbose', 1, 'FitMethod', 'REML'); AIC_values = [omnibus_model.mdlC1.ModelCriterion.AIC omnibus_model.mdlC2.ModelCriterion.AIC ... omnibus_model.mdlC3.ModelCriterion.AIC]; mdlC_AIC = find(AIC_values == min(AIC_values)); % check if the optimal mdlB is better than the optimal mdlA % mdlC_AIC - mdlB_AIC % C2 wins % identify winner of omnibus modelling omnibus_model.winner = omnibus_model.mdlC2 anova(omnibus_model.winner, 'df', 's') %%% analysis of isolated target data targ_tbl = tbl(tbl.condition(:,1) == 't',:); targ_model.base = fitlme(targ_tbl, 'N2pc ~ RT + (RT|subjects) ', 'Verbose', 1, 'FitMethod', 'REML'); targ_model.mdlA1 = fitlme(targ_tbl, ['N2pc ~ RT + (RT|subjects) + ' ... '(RT|targ_cat:targ_index) + (RT|dist_cat:dist_index)'], 'Verbose', 1, 'FitMethod', 'REML'); targ_model.mdlA2 = fitlme(targ_tbl, ['N2pc ~ RT + (RT|subjects) + ' ... '(RT|targ_cat:targ_index) + (1|dist_cat:dist_index)'], 'Verbose', 1, 'FitMethod', 'REML'); % A2 wins targ_model.mdlB1 = fitlme(targ_tbl, ['N2pc ~ RT + (RT|subjects) + ' ... '(1|targ_cat:targ_index) + (1|dist_cat:dist_index)'], 'Verbose', 1, 'FitMethod', 'REML'); % B1 wins targ_model.mdlC1 = fitlme(targ_tbl, ['N2pc ~ RT + (1|subjects) + (1|targ_cat:targ_index) + ' ... '(1|dist_cat:dist_index)'], 'Verbose', 1, 'FitMethod', 'REML'); % C1 wins targ_model.winner = targ_model.mdlC1; anova(targ_model.winner, 'df', 's') %%% analysis of isolated distractor data dist_tbl = tbl(tbl.condition(:,1) == 'd',:); dist_model.base = fitlme(dist_tbl, 'N2pc ~ RT + (RT|subjects) ', 'Verbose', 1, 'FitMethod', 'REML'); dist_model.mdlA1 = fitlme(dist_tbl, ['N2pc ~ RT + (RT|subjects) ' ... '+ (RT|targ_cat:targ_index) + (RT|dist_cat:dist_index)'], 'Verbose', 1, 'FitMethod', 'REML'); dist_model.mdlA2 = fitlme(dist_tbl, ['N2pc ~ RT + (RT|subjects) ' ... '+ (RT|targ_cat:targ_index) + (1|dist_cat:dist_index)'], 'Verbose', 1, 'FitMethod', 'REML'); % A2 wins dist_model.mdlB1 = fitlme(dist_tbl, ['N2pc ~ RT + (RT|subjects) + (1|targ_cat:targ_index) + ' ... '(1|dist_cat:dist_index)'], 'Verbose', 1, 'FitMethod', 'REML'); % B1 wins dist_model.mdlC1 = fitlme(dist_tbl, ['N2pc ~ RT + (1|subjects) + (1|targ_cat:targ_index) + ' ... '(1|dist_cat:dist_index)'], 'Verbose', 1, 'FitMethod', 'REML') % C1 wins dist_model.winner = dist_model.mdlC1; anova(dist_model.winner, 'df', 's') save('models.mat', 'omnibus_model', 'targ_model', 'dist_model'); end targ_tbl = tbl(tbl.condition(:,1) == 't',:); dist_tbl = tbl(tbl.condition(:,1) == 'd',:); %%% eta squared measure of effect size for RT in target model % Correll, Mellinger, & Pedersen, 2021, Behavior Research Methods % Snijders & Bosker, 2011 (SBX) method % Note this method is problematic if random effects are more complicated % than intercepts; see Correll et al. reference. target = targ_model.winner; alternative = fitlme(targ_tbl, 'N2pc ~ 1 + (1|subjects) + (1|targ_cat:targ_index) + (1|dist_cat:dist_index) ', 'Verbose', 1, 'FitMethod', 'REML') eta_squared_targ_model = (alternative.SSE - target.SSE) ./ alternative.SSE; % alternate calculation % target_est = fitted(target); % target_res = targ_tbl.N2pc - target_est; % alternative_est = fitted(alternative); % alternative_res = targ_tbl.N2pc - alternative_est; % diff_res = var(alternative_res) - var(target_res); % eta_squared = diff_res ./ var(alternative_res); %%% check assumptions % linearity - is checked in plot of predictor variable against outcome % variable tbl.estimated_N2pc = fitted(omnibus_model.winner); tbl.residual_N2pc = tbl.N2pc - tbl.estimated_N2pc; % example subject temp_tbl1 = tbl(tbl.subjects(:,1) == 6 & tbl.condition(:,1) == 't',:); temp_tbl2 = tbl(tbl.subjects(:,1) == 6 & tbl.condition(:,1) == 'd',:); figure; hold on; scatter(temp_tbl1.RT, temp_tbl1.estimated_N2pc, 'b.') lsline; scatter(temp_tbl2.RT, temp_tbl2.estimated_N2pc, 'r.') lsline; xlim([0 2200]) title('s6 linearity') % all subjects temp_tbl1 = tbl(tbl.condition(:,1) == 't',:); temp_tbl2 = tbl(tbl.condition(:,1) == 'd',:); figure; hold on; scatter(temp_tbl1.RT, temp_tbl1.estimated_N2pc, 'b.') lsline; scatter(temp_tbl2.RT, temp_tbl2.estimated_N2pc, 'r.') lsline; xlim([0 2200]) title('all linearity') %%% heteroscedasticity - is checked in plot of residuals to data fit, % residuals to predictor figure; plotResiduals(omnibus_model.winner, 'fitted') figure; hold on; scatterhist(tbl.residual_N2pc, tbl.estimated_N2pc, 'Group', tbl.condition, 'kernel', 'on', 'Marker', '.') title('normality of residuals') xlabel('residuals') ylabel('N2pc estimates') figure; scatterhist( tbl.RT, tbl.residual_N2pc, 'Group', tbl.condition, 'kernel', 'on', 'Marker', '.') title('N2pc residual vs RT predictor') % there's a teardrop shape, but this is a product of change in density, not % variance [upper_sd, lower_sd] = deal([]); for sub_num = 1:44 temp = tbl(tbl.subjects == sub_num,:); upper_sd(sub_num) = var(tbl.residual_N2pc(temp.RT > median(temp.RT))) lower_sd(sub_num) = var(tbl.residual_N2pc(temp.RT <= median(temp.RT))) end mean(upper_sd) mean(lower_sd) %%% normality of residuals - is checked in PDFs above, in QQ plot figure; hqqplot = qqplot(tbl.residual_N2pc) set(hqqplot, 'marker', '.') %%% normality of random effects [B Bnames] = randomEffects(omnibus_model.winner) figure; subject_intercept = B(1:2:88); subplot(4,3,1) hist(subject_intercept, 20); xlim([-3 3]); ylim([0 10]); title('subject intercept'); subject_condition = B(2:2:88); subplot(4,3,4) hist(subject_condition, 20); xlim([-3 3]); title('subject condition'); targindex_intercept = B(89:2:408); subplot(4,3,2) hist(targindex_intercept, 20); xlim([-3 3]); title('target intercept'); targindex_condition = B(90:2:408); subplot(4,3,5) hist(targindex_condition, 20); xlim([-3 3]); title('target condition') distindex_intercept = B(409:2:end); subplot(4,3,3) hist(distindex_intercept, 20); xlim([-3 3]); title('dist intercept'); distindex_condition = B(410:2:end); subplot(4,3,6) hist(distindex_condition, 20); xlim([-3 3]); title('dist condition') [B Bnames] = randomEffects(targ_model.winner) subject_intercept = B(1:44); subplot(4,3,7) hist(subject_intercept, 20); xlim([-3 3]); ylim([0 10]); title('target model - subject intercept'); targindex_intercept = B(45:200); subplot(4,3,8) hist(targindex_intercept, 20); xlim([-3 3]); title('target model - targ intercept'); distindex_intercept = B(201:end); subplot(4,3,9) hist(distindex_intercept, 20); xlim([-3 3]); title('target model - dist intercept'); [B Bnames] = randomEffects(dist_model.winner) subject_intercept = B(1:44); subplot(4,3,10) hist(subject_intercept, 20); xlim([-3 3]); ylim([0 10]); title('dist model - subject intercept'); targindex_intercept = B(45:200); subplot(4,3,11) hist(targindex_intercept, 20); xlim([-3 3]); title('dist model - targ intercept'); distindex_intercept = B(201:end); subplot(4,3,12) hist(distindex_intercept, 20); xlim([-3 3]); title('dist model - dist intercept'); targ_ML = fitlme(targ_tbl, 'N2pc ~ RT + (1|subjects) + (1|targ_cat:targ_index) + (1|dist_cat:dist_index) ', 'Verbose', 1, 'FitMethod', 'ML') targ_ML_nodist = fitlme(targ_tbl, 'N2pc ~ RT + (1|subjects) + (1|targ_cat:targ_index) ', 'Verbose', 1, 'FitMethod', 'ML') targ_ML_notarg = fitlme(targ_tbl, 'N2pc ~ RT + (1|subjects) + (1|dist_cat:dist_index) ', 'Verbose', 1, 'FitMethod', 'ML') compare(targ_ML_notarg, targ_ML)