Compare Models ********** Now we're in the home stretch. We've gone through all of these steps to get here: the point where we can now draw conclusions about human decision-making and how it is influenced by their social preferences. Lesson ================ .. article-info:: :avatar: UCLA_Suit.jpg :avatar-link: https://www.decisionneurosciencelab.com/elijahgalvan :author: Elijah Galvan :date: September 1, 2023 :read-time: 8 min read :class-container: sd-p-2 sd-outline-muted sd-rounded-1 Goal During this Stage --------------- We're going to compare models against each other now. While we've already determined what the best model is, we haven't yet tested the hypotheses that it embodies. How to Achieve this Goal ------------ .. dropdown:: Hypothesis Testing via Model Comparison .. tab-set:: .. tab-item:: Plain English We're going to compare MFIs for various models within subjects using uncorrected paired t-tests. Remember that we have to exclude any subjects for whom the model perfectly explains behavior. Also, if we hypothesized that our favored model would outperform other models, we should use directional tests. .. Note:: You should never correct for multiple comparisons when testing computational models against each other. Your model fits the data best: you're testing the best model, and the best model only, against all of the other models in your set. One of those models in your model set will have the lowest lower boundary of their confidence interval - this is the model we want to test against (this is usually the second-best performing model but not necessarily). There is no reason to compromise our power when we are making multiple comparisons only to facilitate testing your model against the best challenger. .. tab-item:: R :: excluded1 = which(subjectData$modelSS == 0 | altSubjectData$modelAlt1SS == 0) excluded2 = which(subjectData$modelSS == 0 | altSubjectData$mode2Alt1SS == 0) t.test(subjectData$modelAIC[-excluded1], altSubjectData$modelAlt1AIC[-excluded1], paired = T, alternative = 'less') #favored model should be less than the other model (i.e. better model fit) t.test(subjectData$modelAIC[-excluded2], altSubjectData$modelAlt2AIC[-excluded2], paired = T, alternative = 'less') .. tab-item:: MatLab :: excluded1 = find(subjectData.modelSS == 0 | altSubjectData.modelAlt1SS == 0); excluded2 = find(subjectData.modelSS == 0 | altSubjectData.modelAlt2SS == 0); [~, p1] = ttest(subjectData.modelAIC(setdiff(1:end, excluded1)), altSubjectData.modelAlt1AIC(setdiff(1:end, excluded1)), 'Tail', 'left'); [~, p2] = ttest(subjectData.modelAIC(setdiff(1:end, excluded2)), altSubjectData.modelAlt2AIC(setdiff(1:end, excluded2)), 'Tail', 'left'); .. tab-item:: Python :: from scipy.stats import ttest_rel excluded1 = np.where((subjectData['modelSS'] == 0) | (altSubjectData['modelAlt1SS'] == 0))[0] excluded2 = np.where((subjectData['modelSS'] == 0) | (altSubjectData['modelAlt2SS'] == 0))[0] _, p1 = ttest_rel(subjectData['modelAIC'][np.setdiff1d(np.arange(len(subjectData)), excluded1)], altSubjectData['modelAlt1AIC'][np.setdiff1d(np.arange(len(altSubjectData)), excluded1)], alternative='less') _, p2 = ttest_rel(subjectData['modelAIC'][np.setdiff1d(np.arange(len(subjectData)), excluded2)], altSubjectData['modelAlt2AIC'][np.setdiff1d(np.arange(len(altSubjectData)), excluded2)], alternative='less') Tutorials ========== Tutorial 1 - van Baar, Chang, & Sanfey, 2019 ---------------------- .. dropdown:: Hypothesis Testing via Model Comparison .. tab-set:: .. tab-item:: R :: t.test(subjectData$modelAIC, altSubjectData$modelAICGreed, paired = T) #negative mean difference/t value means that first term is lower than second term and lower MFI is better t.test(subjectData$modelAIC, altSubjectData$modelAICGuilt, paired = T) t.test(subjectData$modelAIC, altSubjectData$modelAICInequity, paired = T) library(ggsignif) aic = c(mean(subjectData$modelAIC), mean(altSubjectData$modelAICGreed), mean(altSubjectData$modelAICGuilt), mean(altSubjectData$modelAICInequity)) qplot(y = aic, x = as.factor(c('Moral Strategies Model', 'Greed Model', 'Guilt Model', 'Inequity Model')), fill = as.factor(c('Moral Strategies Model', 'Greed Model', 'Guilt Model', 'Inequity Model')), color = '', geom = 'col') + labs(x = 'Model', y = 'AIC', fill = NULL, color = NULL) + theme_minimal() + scale_color_manual(values = c(rgb(0, 0, 0, maxColorValue = 255))) + scale_fill_manual(values = c(rgb(132.5, 132.5, 132.5, maxColorValue = 255), rgb(132.5, 132.5, 132.5, maxColorValue = 255), rgb(132.5, 132.5, 132.5, maxColorValue = 255), rgb(218, 165, 32, maxColorValue = 255))) + geom_signif(comparisons = list(c('Moral Strategies Model', 'Greed Model')), y = 341, textsize = 0)+ geom_signif(comparisons = list(c('Moral Strategies Model', 'Guilt Model')), y = 328, textsize = 0)+ geom_signif(comparisons = list(c('Moral Strategies Model', 'Inequity Model')), y = 315, textsize = 0) + annotate("text", x = 2.5, label = "***", y = 358, size = 4)+ # 3 for p < 0.001, 2 for p < 0.01, 1 for p < 0.05 usually annotate("text", x = 3, label = "***", y = 345, size = 4)+ annotate("text", x = 3.5, label = "***", y = 332, size = 4) .. tab-item:: MatLab :: ttest(subjectData.modelAIC, altSubjectData.modelAICGreed, 'Paired', true); ttest(subjectData.modelAIC, altSubjectData.modelAICGuilt, 'Paired', true); ttest(subjectData.modelAIC, altSubjectData.modelAICInequity, 'Paired', true); aic = [mean(subjectData.modelAIC), mean(altSubjectData.modelAICGreed), mean(altSubjectData.modelAICGuilt), mean(altSubjectData.modelAICInequity)]; figure; bar(categorical({'Moral Strategies Model', 'Greed Model', 'Guilt Model', 'Inequity Model'}), aic, 'FaceColor', [0.85 0.85 0.85], 'EdgeColor', 'none'); ylabel('AIC'); xlabel('Model'); ylim([300 360]); hold on; text(2.5, 358, '***', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 14); text(3, 345, '***', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 14); text(3.5, 332, '***', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 14); hold off; .. tab-item:: Python :: ttest_rel(subjectData['modelAIC'], altSubjectData['modelAICGreed']) ttest_rel(subjectData['modelAIC'], altSubjectData['modelAICGuilt']) ttest_rel(subjectData['modelAIC'], altSubjectData['modelAICInequity']) aic = [np.mean(subjectData['modelAIC']), np.mean(altSubjectData['modelAICGreed']), np.mean(altSubjectData['modelAICGuilt']), np.mean(altSubjectData['modelAICInequity'])] plt.figure(); plt.bar(['Moral Strategies Model', 'Greed Model', 'Guilt Model', 'Inequity Model'], aic, color=[0.85, 0.85, 0.85], edgecolor='none'); plt.ylabel('AIC'); plt.xlabel('Model'); plt.ylim([300, 360]); plt.text(2.5, 358, '***', ha='center', va='bottom', fontsize=14); plt.text(3, 345, '***', ha='center', va='bottom', fontsize=14); plt.text(3.5, 332, '***', ha='center', va='bottom', fontsize=14); plt.show(); Tutorial 2 - Galvan & Sanfey, 2024 ------------------- .. dropdown:: Hypothesis Testing via Model Comparison .. tab-set:: .. tab-item:: R :: alpha = 0.05/length(conditions) # Merit Condition t.test(subjectData$AICMerit[-excludedM], altSubjectData$AICPMerit[-excludedM], paired = T, conf.level = (1 - alpha), alternative = 'less') t.test(subjectData$AICMerit[-excludedM], altSubjectData$AICLMerit[-excludedM], paired = T, conf.level = (1 - alpha), alternative = 'less') t.test(subjectData$AICMerit[-excludedM], altSubjectData$AICTMerit[-excludedM], paired = T, conf.level = (1 - alpha), alternative = 'less') t.test(subjectData$AICMerit[-excludedM], altSubjectData$AICPLMerit[-excludedM], paired = T, conf.level = (1 - alpha), alternative = 'less') t.test(subjectData$AICMerit[-excludedM], altSubjectData$AICPTMerit[-excludedM], paired = T, conf.level = (1 - alpha), alternative = 'less') t.test(subjectData$AICMerit[-excludedM], altSubjectData$AICLTMerit[-excludedM], paired = T, conf.level = (1 - alpha), alternative = 'less') # Corruption Condition t.test(subjectData$AICEntitlement[-excludedE], altSubjectData$AICPEntitlement[-excludedE], paired = T, conf.level = (1 - alpha), alternative = 'less') t.test(subjectData$AICEntitlement[-excludedE], altSubjectData$AICLEntitlement[-excludedE], paired = T, conf.level = (1 - alpha), alternative = 'less') t.test(subjectData$AICEntitlement[-excludedE], altSubjectData$AICTEntitlement[-excludedE], paired = T, conf.level = (1 - alpha), alternative = 'less') t.test(subjectData$AICEntitlement[-excludedE], altSubjectData$AICPLEntitlement[-excludedE], paired = T, conf.level = (1 - alpha), alternative = 'less') t.test(subjectData$AICEntitlement[-excludedE], altSubjectData$AICPTEntitlement[-excludedE], paired = T, conf.level = (1 - alpha), alternative = 'less') t.test(subjectData$AICEntitlement[-excludedE], altSubjectData$AICLTEntitlement[-excludedE], paired = T, conf.level = (1 - alpha), alternative = 'less') # Corruption Condition t.test(subjectData$AICCorruption[-excludedC], altSubjectData$AICPCorruption[-excludedC], paired = T, conf.level = (1 - alpha), alternative = 'less') t.test(subjectData$AICCorruption[-excludedC], altSubjectData$AICLCorruption[-excludedC], paired = T, conf.level = (1 - alpha), alternative = 'less') t.test(subjectData$AICCorruption[-excludedC], altSubjectData$AICTCorruption[-excludedC], paired = T, conf.level = (1 - alpha), alternative = 'less') t.test(subjectData$AICCorruption[-excludedC], altSubjectData$AICPLCorruption[-excludedC], paired = T, conf.level = (1 - alpha), alternative = 'less') t.test(subjectData$AICCorruption[-excludedC], altSubjectData$AICPTCorruption[-excludedC], paired = T, conf.level = (1 - alpha), alternative = 'less') t.test(subjectData$AICCorruption[-excludedC], altSubjectData$AICLTCorruption[-excludedC], paired = T, conf.level = (1 - alpha), alternative = 'less') # Luck Condition t.test(subjectData$AICLuck[-excludedL], altSubjectData$AICPLuck[-excludedL], paired = T, conf.level = (1 - alpha), alternative = 'less') t.test(subjectData$AICLuck[-excludedL], altSubjectData$AICLLuck[-excludedL], paired = T, conf.level = (1 - alpha), alternative = 'less') t.test(subjectData$AICLuck[-excludedL], altSubjectData$AICTLuck[-excludedL], paired = T, conf.level = (1 - alpha), alternative = 'less') t.test(subjectData$AICLuck[-excludedL], altSubjectData$AICPLLuck[-excludedL], paired = T, conf.level = (1 - alpha), alternative = 'less') t.test(subjectData$AICLuck[-excludedL], altSubjectData$AICPTLuck[-excludedL], paired = T, conf.level = (1 - alpha), alternative = 'less') t.test(subjectData$AICLuck[-excludedL], altSubjectData$AICLTLuck[-excludedL], paired = T, conf.level = (1 - alpha), alternative = 'less') .. tab-item:: MatLab :: alpha = 0.05 / length(conditions); % Merit Condition for i = 1:length(conditions) ttest(subjectData.AICMerit(~excludedM), altSubjectData.AICPMerit(~excludedM), 'Tail', 'left', 'Alpha', alpha); ttest(subjectData.AICMerit(~excludedM), altSubjectData.AICLMerit(~excludedM), 'Tail', 'left', 'Alpha', alpha); ttest(subjectData.AICMerit(~excludedM), altSubjectData.AICTMerit(~excludedM), 'Tail', 'left', 'Alpha', alpha); ttest(subjectData.AICMerit(~excludedM), altSubjectData.AICPLMerit(~excludedM), 'Tail', 'left', 'Alpha', alpha); ttest(subjectData.AICMerit(~excludedM), altSubjectData.AICPTMerit(~excludedM), 'Tail', 'left', 'Alpha', alpha); ttest(subjectData.AICMerit(~excludedM), altSubjectData.AICLTMerit(~excludedM), 'Tail', 'left', 'Alpha', alpha); end % Entitlement Condition for i = 1:length(conditions) ttest(subjectData.AICEntitlement(~excludedE), altSubjectData.AICPEntitlement(~excludedE), 'Tail', 'left', 'Alpha', alpha); ttest(subjectData.AICEntitlement(~excludedE), altSubjectData.AICLEntitlement(~excludedE), 'Tail', 'left', 'Alpha', alpha); ttest(subjectData.AICEntitlement(~excludedE), altSubjectData.AICTEntitlement(~excludedE), 'Tail', 'left', 'Alpha', alpha); ttest(subjectData.AICEntitlement(~excludedE), altSubjectData.AICPLEntitlement(~excludedE), 'Tail', 'left', 'Alpha', alpha); ttest(subjectData.AICEntitlement(~excludedE), altSubjectData.AICPTEntitlement(~excludedE), 'Tail', 'left', 'Alpha', alpha); ttest(subjectData.AICEntitlement(~excludedE), altSubjectData.AICLTEntitlement(~excludedE), 'Tail', 'left', 'Alpha', alpha); end % Corruption Condition for i = 1:length(conditions) ttest(subjectData.AICCorruption(~excludedC), altSubjectData.AICPCorruption(~excludedC), 'Tail', 'left', 'Alpha', alpha); ttest(subjectData.AICCorruption(~excludedC), altSubjectData.AICLCorruption(~excludedC), 'Tail', 'left', 'Alpha', alpha); ttest(subjectData.AICCorruption(~excludedC), altSubjectData.AICTCorruption(~excludedC), 'Tail', 'left', 'Alpha', alpha); ttest(subjectData.AICCorruption(~excludedC), altSubjectData.AICPLCorruption(~excludedC), 'Tail', 'left', 'Alpha', alpha); ttest(subjectData.AICCorruption(~excludedC), altSubjectData.AICPTCorruption(~excludedC), 'Tail', 'left', 'Alpha', alpha); ttest(subjectData.AICCorruption(~excludedC), altSubjectData.AICLTCorruption(~excludedC), 'Tail', 'left', 'Alpha', alpha); end % Luck Condition for i = 1:length(conditions) ttest(subjectData.AICLuck(~excludedL), altSubjectData.AICPLuck(~excludedL), 'Tail', 'left', 'Alpha', alpha); ttest(subjectData.AICLuck(~excludedL), altSubjectData.AICLLuck(~excludedL), 'Tail', 'left', 'Alpha', alpha); ttest(subjectData.AICLuck(~excludedL), altSubjectData.AICTLuck(~excludedL), 'Tail', 'left', 'Alpha', alpha); ttest(subjectData.AICLuck(~excludedL), altSubjectData.AICPLLuck(~excludedL), 'Tail', 'left', 'Alpha', alpha); ttest(subjectData.AICLuck(~excludedL), altSubjectData.AICPTLuck(~excludedL), 'Tail', 'left', 'Alpha', alpha); ttest(subjectData.AICLuck(~excludedL), altSubjectData.AICLTLuck(~excludedL), 'Tail', 'left', 'Alpha', alpha); end .. tab-item:: Python :: alpha <- 0.05 / length(conditions) # Merit Condition t.test(subjectData$AICMerit[-excludedM], altSubjectData$AICPMerit[-excludedM], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') t.test(subjectData$AICMerit[-excludedM], altSubjectData$AICLMerit[-excludedM], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') t.test(subjectData$AICMerit[-excludedM], altSubjectData$AICTMerit[-excludedM], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') t.test(subjectData$AICMerit[-excludedM], altSubjectData$AICPLMerit[-excludedM], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') t.test(subjectData$AICMerit[-excludedM], altSubjectData$AICPTMerit[-excludedM], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') t.test(subjectData$AICMerit[-excludedM], altSubjectData$AICLTMerit[-excludedM], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') # Entitlement Condition t.test(subjectData$AICEntitlement[-excludedE], altSubjectData$AICPEntitlement[-excludedE], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') t.test(subjectData$AICEntitlement[-excludedE], altSubjectData$AICLEntitlement[-excludedE], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') t.test(subjectData$AICEntitlement[-excludedE], altSubjectData$AICTEntitlement[-excludedE], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') t.test(subjectData$AICEntitlement[-excludedE], altSubjectData$AICPLEntitlement[-excludedE], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') t.test(subjectData$AICEntitlement[-excludedE], altSubjectData$AICPTEntitlement[-excludedE], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') t.test(subjectData$AICEntitlement[-excludedE], altSubjectData$AICLTEntitlement[-excludedE], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') # Corruption Condition t.test(subjectData$AICCorruption[-excludedC], altSubjectData$AICPCorruption[-excludedC], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') t.test(subjectData$AICCorruption[-excludedC], altSubjectData$AICLCorruption[-excludedC], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') t.test(subjectData$AICCorruption[-excludedC], altSubjectData$AICTCorruption[-excludedC], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') t.test(subjectData$AICCorruption[-excludedC], altSubjectData$AICPLCorruption[-excludedC], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') t.test(subjectData$AICCorruption[-excludedC], altSubjectData$AICPTCorruption[-excludedC], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') t.test(subjectData$AICCorruption[-excludedC], altSubjectData$AICLTCorruption[-excludedC], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') # Luck Condition t.test(subjectData$AICLuck[-excludedL], altSubjectData$AICPLuck[-excludedL], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') t.test(subjectData$AICLuck[-excludedL], altSubjectData$AICLLuck[-excludedL], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') t.test(subjectData$AICLuck[-excludedL], altSubjectData$AICTLuck[-excludedL], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') t.test(subjectData$AICLuck[-excludedL], altSubjectData$AICPLLuck[-excludedL], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') t.test(subjectData$AICLuck[-excludedL], altSubjectData$AICPTLuck[-excludedL], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') t.test(subjectData$AICLuck[-excludedL], altSubjectData$AICLTLuck[-excludedL], paired = TRUE, conf.level = 1 - alpha, alternative = 'less') Tutorial 3 - Crockett et al., 2014 ------------------- .. dropdown:: Hypothesis Testing via Model Comparison .. tab-set:: .. tab-item:: R :: t.test(subjectData$BIC[-20], altSubjectData$BICM2[-20], paired = T) t.test(subjectData$BIC[-20], altSubjectData$BICM3[-20], paired = T) t.test(subjectData$BIC[-20], altSubjectData$BICM4[-20], paired = T) t.test(subjectData$BIC[-20], altSubjectData$BICM5[-20], paired = T) .. tab-item:: MatLab :: .. tab-item:: Python :: Tutorial 4 - Li et al., 2022 ------------------- .. dropdown:: Hypothesis Testing via Model Comparison .. tab-set:: .. tab-item:: R :: t.test(altSubjectData$BIC_M5, subjectData$BIC, paired = T) t.test(altSubjectData$BIC_M5, altSubjectData$BIC_M1, paired = T) t.test(altSubjectData$BIC_M5, altSubjectData$BIC_M2, paired = T) t.test(altSubjectData$BIC_M5, altSubjectData$BIC_M3, paired = T) t.test(altSubjectData$BIC_M5, altSubjectData$BIC_M4, paired = T) t.test(altSubjectData$BIC_M5, altSubjectData$BIC_M6, paired = T) .. tab-item:: MatLab :: .. tab-item:: Python ::