Identify the Best Model
Lesson
Elijah Galvan
September 1, 2023
15 min read
Goal During this Stage
To determine which model is best.
Note
Now if you’ve been following along, you’ve probably realized that we’ve only talked so far about a single model - I’ll briefly explain why. We use Social Utility Models to encapsulate individual differences in value based preferences and exisiting research on this topic essentially shows that people are usually different in most conceivable ways. Thus, we always try to design experiments around trying to see this individual differences in behavior and our “favored” model - the one we think will best fit the data - is the one which captures all of these differences. Since the favored model must capture all meaningful potential differences, it is therefore the most complex model that we will test.
So why are we reiterating this now? Well, the answer is that any other types of models you test on value-based decision-making behavior in non-forced, non-probabilistic, social choice tasks are almost always going to perform much worse. These tasks are designed to elicit consistent, well-defined preferences with a very high signal to noise ratio: therefore, models capturing other aspects of the decision-making process other than these preferences are theoretically and pragmatically inferior in these scenarios. All of this to say, we are going to test this favored model - which is the most complex model - against models which are simplifications of this model. Thus, we can use the same Construct Value computations, removing one of these Construct Value computations from the utility equation and a Free Parameter per model.
How to Achieve this Goal
Create Utility Equations for our Alternative Models
Generally, you should first enumerate all of the reasonable potential alternative models. This means that if you have a model with 3 different Construct Value terms (as most models that we will talk about indeed do) then the alternative model set would include models with only 2 or 1 Construct Value terms. Are any models not worth testing? Then don’t include them. For every model we think is worth testing, let’s talk about how you’re going to go about testing them.
Any models that we want to test with one Construct Value term are really, really simple - Utility = Construct Value. Notice that there are no Free Parameters in this Equation - since only one Construct Value determines Utility, there is no way that it can capture individual differences.
The models that we want to test with two Construct Value terms are also pretty simple - we want to capture how people prioritize one norm over another. Thus, we express this as Utility = (Construct Value 1 x Free Parameter) + (Construct Value 2 x (1 - Free Parameter)).
Preallocating and Defining Functions
We’re going to need several items for this section:
Utility Functions for each each Model
Objective Functions for each Utility Equation with Free Parameters
New Optimizer Inputs for each Utility Equation with Free Parameters
Outputs for the Alternative Models
utility_alt1 = function(construct1){
return(utility)
}
utility_alt2 = function(construct1, construct2, parameter1){
return(utility)
}
#alternative model with 1 construct has 0 free parameters: doesn't need an objective function
obj_function_alt2 = function(param, decisions, method = "OLS") {
Parameter1 = param[1]
predicted_utility = vector('numeric', length(trialList[,1]))
chosen = decisions + 1
for (k in 1:length(trialList[,1])){
IV = trialList[k, 1]
Constant = trialList[k, 2]
Choices = seq(0, (I * M), 1)
Utility = vector('numeric', length(Choices))
for (n in 1:length(Choices)){
Utility[n] = utility_alt2(Parameter1, construct1(IV, Constant, Choices[n]))
}
predicted_utility[k] = max(Utility)
observed_utility[k] = Utility[chosen[k]]
}
if (method == "OLS"){
return(sum((predicted_utility - observed_utility)**2))
} else if (method == "MLE"){
return(-1 * sum(dnorm(observed_utility, mean = predicted_utility, sd = sd, log = TRUE)))
}
}
initial_param_alt2 = #something
lower_bound_alt2 = #something
upper_bound_alt2 = #something
altSubjectData = data.frame()
function utility = utility_alt1(construct1)
% Your utility calculation for utility_alt1 goes here
utility = 0; % Placeholder value, replace with actual calculation
end
function utility = utility_alt2(construct1, construct2, parameter1)
% Your utility calculation for utility_alt2 goes here
utility = 0; % Placeholder value, replace with actual calculation
end
function obj_value = obj_function_alt2(param, decisions, method)
Parameter1 = param(1);
trialList = []; % Define trialList here
predicted_utility = zeros(size(trialList, 1), 1);
chosen = decisions + 1;
for k = 1:size(trialList, 1)
IV = trialList(k, 1);
Constant = trialList(k, 2);
Choices = 0:(I * M); % Define I and M
Utility = zeros(size(Choices));
for n = 1:length(Choices)
Utility(n) = utility_alt2(Parameter1, construct1(IV, Constant, Choices(n)));
end
predicted_utility(k) = max(Utility);
observed_utility(k) = Utility(chosen(k));
end
if strcmp(method, 'OLS')
obj_value = sum((predicted_utility - observed_utility).^2);
elseif strcmp(method, 'MLE')
obj_value = -sum(log(normpdf(observed_utility, predicted_utility, sd)));
end
end
% Define initial_param_alt2, lower_bound_alt2, and upper_bound_alt2 here
altSubjectData = table(); % Create an empty table for altSubjectData
from scipy.stats import norm
def utility_alt1(construct1):
# Your utility calculation for utility_alt1 goes here
utility = 0 # Placeholder value, replace with actual calculation
return utility
def utility_alt2(construct1, construct2, parameter1):
# Your utility calculation for utility_alt2 goes here
utility = 0 # Placeholder value, replace with actual calculation
return utility
def obj_function_alt2(param, decisions, method="OLS"):
Parameter1 = param[0] # MATLAB-style indexing
trialList = np.array([]) # Define trialList here
predicted_utility = np.zeros(len(trialList))
chosen = decisions + 1
for k in range(len(trialList)):
IV = trialList[k, 0]
Constant = trialList[k, 1]
Choices = np.arange(0, I * M + 1) # Define I and M
Utility = np.zeros(len(Choices))
for n in range(len(Choices)):
Utility[n] = utility_alt2(Parameter1, construct1(IV, Constant, Choices[n]))
predicted_utility[k] = np.max(Utility)
observed_utility[k] = Utility[chosen[k]]
if method == "OLS":
obj_value = np.sum((predicted_utility - observed_utility)**2)
elif method == "MLE":
obj_value = -np.sum(np.log(norm.pdf(observed_utility, predicted_utility, sd)))
return obj_value
# Define initial_param_alt2, lower_bound_alt2, and upper_bound_alt2 here
altSubjectData = pd.DataFrame() # Create an empty DataFrame for altSubjectData
Recover Free Parameters for Each Alternative Model, Per Subject
For the alternative models with Free Parameters, we’ll need to recover these Free Parameters in order to generate model predictions. Let’s do that quickly in the same way that we did for the other model, leaving a demand to subsequently determine model predictions.
Note
Models with only one Free Parameter may require a different optimzer than models with multiple Free Parameters. However, you won’t have to change anything about how your objective function works or anything like that so don’t worry!
for (i in 1:length(included_subjects)){
datafile = paste(parentfolder, included_subjects[i], restoffilepath, sep = '') # produces a character vector 'parentfolder/included_subjects[i]**.filetype'
df = read.csv2(datafile)
reorder = df$trialsTask.thisIndex + 1
result_alt2 = optim(obj_function_alt2, par = initial_param_alt2, lower = lower_bound_alt2, upper = upper_bound_alt2, decisions = df$Decisions)
# Determine Predictions
}
for i = 1:length(included_subjects)
datafile = strcat(parentfolder, included_subjects{i}, restoffilepath); % produces a character vector 'parentfolder/included_subjects{i}**.filetype'
df = readtable(datafile);
reorder = df.trialsTask.thisIndex + 1;
result_alt2 = fmincon(@obj_function_alt2, initial_param_alt2, [], [], [], [], lower_bound_alt2, upper_bound_alt2, [], optimset('Display', 'off'), df.Decisions);
% Determine Predictions
end
for i in range(len(included_subjects)):
datafile = parentfolder + included_subjects[i] + restoffilepath # produces a character vector 'parentfolder/included_subjects[i]**.filetype'
df = pd.read_csv(datafile, sep='\t')
reorder = df['trialsTask.thisIndex'] + 1
result_alt2 = optimize.minimize(obj_function_alt2, initial_param_alt2, bounds=list(zip(lower_bound_alt2, upper_bound_alt2)), args=(df['Decisions'],))
# Determine Predictions
Determine Predicted Decisions for Each Alternative Model, Per Subject
Now, we are going to answer the Determine Predictions demand placed on us. We have found the Subject’s Free Parameters so we need to specifically know what it is that our model predicts that they will do. In the previous step, we could have cut a corner and gotten the predictions from the closest point we simulated data for. In all likelihood, the model predictions would be indistinguishable from these, but for the sake of being punctual let’s get these predictions in the same way we did with our favored model!
for (i in 1:length(included_subjects)){
datafile = paste(parentfolder, included_subjects[i], restoffilepath, sep = '') # produces a character vector 'parentfolder/included_subjects[i]**.filetype'
df = read.csv2(datafile)
reorder = df$trialsTask.thisIndex + 1
result_alt2 = optim(obj_function_alt2, par = initial_param_alt2, lower = lower_bound_alt2, upper = upper_bound_alt2, decisions = df$Decisions)
#Just Added
df$PredictionAlt1 = vector('numeric')
df$PredictionAlt2 = df$PredictionAlt1
for (k in 1:length(df$Decisions)){
UtilityAlt1 = vector('numeric', length(Choices))
UtilityAlt2 = vector('numeric', length(Choices))
for (n in 1:length(Choices)){
UtilityAlt1[n] = utility_alt1(construct1 = construct1(df$IV[k], df$Constant[k], Choices[n]))
UtilityAlt2[n] = utility_alt2(parameter1 = result_alt2$par[1],
construct1 = construct1(df$IV[k], df$Constant[k], Choices[n]),
construct2 = construct2(df$IV[k], df$Constant[k], Choices[n]))
}
correct_choice_alt1 = which(UtilityAlt1 == max(UtilityAlt1))
correct_choice_alt2 = which(UtilityAlt2 == max(UtilityAlt2))
if (length(correct_choice) > 1){
correct_choice = correct_choice[sample(correct_choice, 1)]
}
df$PredictionAlt1[k] = Choices[correct_choice_alt1]
df$PredictionAlt2[k] = Choices[correct_choice_alt2]
}
model_NLL_Alt1 = -2 * log(sum(dnorm(df$Decision, mean = df$PredictionAlt1)))
model_SS_Alt1 = sum((df$Decision - df$PredictionAlt1)**2)
model_NLL_Alt2 = -2 * log(sum(dnorm(df$Decision, mean = df$PredictionAlt2)))
model_SS_Alt2 = sum((df$Decision - df$PredictionAlt2)**2)
altSubjectData[i, 1:6] = c(included_subjects[i], result_alt2$par[1], model_NLL_Alt1, model_SS_Alt1, model_NLL_Alt2, model_SS_Alt2)
}
colnames(altSubjectData) = c('SubjectID', 'Parameter1_Alt2', 'alt1_modelNLL', 'alt1_modelSS', 'alt2_modelNLL', 'alt2_modelSS')
for i = 1:length(included_subjects)
datafile = strcat(parentfolder, included_subjects{i}, restoffilepath);
df = readtable(datafile);
reorder = df.trialsTask.thisIndex + 1;
result_alt2 = fmincon(@(x) obj_function_alt2(x, df.Decisions), initial_param_alt2, [], [], [], [], lower_bound_alt2, upper_bound_alt2);
% Just Added
df.PredictionAlt1 = zeros(size(df.Decisions));
df.PredictionAlt2 = df.PredictionAlt1;
for k = 1:length(df.Decisions)
UtilityAlt1 = zeros(size(Choices));
UtilityAlt2 = zeros(size(Choices));
for n = 1:length(Choices)
UtilityAlt1(n) = utility_alt1(construct1(df.IV(k), df.Constant(k), Choices(n)));
UtilityAlt2(n) = utility_alt2(result_alt2(1), construct1(df.IV(k), df.Constant(k), Choices(n)), construct2(df.IV(k), df.Constant(k), Choices(n)));
end
[~, correct_choice_alt1] = max(UtilityAlt1);
[~, correct_choice_alt2] = max(UtilityAlt2);
if length(correct_choice_alt1) > 1
correct_choice_alt1 = correct_choice_alt1(randi(length(correct_choice_alt1)));
end
df.PredictionAlt1(k) = Choices(correct_choice_alt1);
df.PredictionAlt2(k) = Choices(correct_choice_alt2);
end
model_NLL_Alt1 = -2 * sum(log(normpdf(df.Decisions, df.PredictionAlt1)));
model_SS_Alt1 = sum((df.Decisions - df.PredictionAlt1).^2);
model_NLL_Alt2 = -2 * sum(log(normpdf(df.Decisions, df.PredictionAlt2)));
model_SS_Alt2 = sum((df.Decisions - df.PredictionAlt2).^2);
altSubjectData(i, 1:6) = [included_subjects{i}, result_alt2(1), model_NLL_Alt1, model_SS_Alt1, model_NLL_Alt2, model_SS_Alt2];
end
altSubjectData.Properties.VariableNames = {'SubjectID', 'Parameter1_Alt2', 'alt1_modelNLL', 'alt1_modelSS', 'alt2_modelNLL', 'alt2_modelSS'};
for i in range(len(included_subjects)):
datafile = parentfolder + included_subjects[i] + restoffilepath
df = pd.read_csv(datafile, delimiter=',')
reorder = df['trialsTask.thisIndex'] + 1
result_alt2 = minimize(lambda x: obj_function_alt2(x, df['Decisions']), initial_param_alt2, bounds=list(zip(lower_bound_alt2, upper_bound_alt2)))
# Just Added
df['PredictionAlt1'] = np.zeros(len(df['Decisions']))
df['PredictionAlt2'] = df['PredictionAlt1']
for k in range(len(df['Decisions'])):
UtilityAlt1 = np.zeros(len(Choices))
UtilityAlt2 = np.zeros(len(Choices))
for n in range(len(Choices)):
UtilityAlt1[n] = utility_alt1(construct1(df['IV'][k], df['Constant'][k], Choices[n]))
UtilityAlt2[n] = utility_alt2(result_alt2.x[0], construct1(df['IV'][k], df['Constant'][k], Choices[n]), construct2(df['IV'][k], df['Constant'][k], Choices[n]))
correct_choice_alt1 = np.argmax(UtilityAlt1)
correct_choice_alt2 = np.argmax(UtilityAlt2)
if len(correct_choice_alt1) > 1:
correct_choice_alt1 = np.random.choice(correct_choice_alt1, 1)
df['PredictionAlt1'][k] = Choices[correct_choice_alt1]
df['PredictionAlt2'][k] = Choices[correct_choice_alt2]
model_NLL_Alt1 = -2 * np.sum(np.log(np.random.normal(df['Decisions'], df['PredictionAlt1'])))
model_SS_Alt1 = np.sum((df['Decisions'] - df['PredictionAlt1'])**2)
model_NLL_Alt2 = -2 * np.sum(np.log(np.random.normal(df['Decisions'], df['PredictionAlt2'])))
model_SS_Alt2 = np.sum((df['Decisions'] - df['PredictionAlt2'])**2)
altSubjectData[i, 1:6] = [included_subjects[i], result_alt2.x[0], model_NLL_Alt1, model_SS_Alt1, model_NLL_Alt2, model_SS_Alt2]
altSubjectData.columns = ['SubjectID', 'Parameter1_Alt2', 'alt1_modelNLL', 'alt1_modelSS', 'alt2_modelNLL', 'alt2_modelSS']
Compute Model Fit Index for Each Subject, for Each Alternative Model
Now that we have the model error - either the sum of squared errors or the negative log likelihood of the real Subjects Decisions versus the model’s predicted Decisions.
altSubjectData$modelAlt1AIC = N * log(altSubjectData$modelSS_Alt1/N) + 2*0
altSubjectData$modelAlt2AIC = N * log(altSubjectData$modelSS_Alt2/N) + 2*1
altSubjectData.modelAlt1AIC = N * log(altSubjectData.modelSS_Alt1/N) + 2*0;
altSubjectData.modelAlt2AIC = N * log(altSubjectData.modelSS_Alt2/N) + 2*1;
altSubjectData['modelAlt1AIC'] = N * log(altSubjectData['modelSS_Alt1']/N) + 2*0
altSubjectData['modelAlt2AIC'] = N * log(altSubjectData['modelSS_Alt2']/N) + 2*1
Compare Model Performance
Now we simply want to identify which model is best. Thus, we’re going to create a vector with the MFI for each model averaged across all Subjects and select the model with the lowest MFI. This approach tells us which model provides the best average fit for Subjects in our sample. Importantly, if any Subjects data are fully explained by the model (i.e. observed Decisions always equal Decisions predicted by the model) then these Subjects must be excluded from your analysis since their MFIs are negative infinity.
Another approach would be to compute MFIs for the entire dataset - this approach does not require that you exclude Subjects from analysis. It is more appropriate to use the first approach if you are focused on individual differences (i.e. trying to characterize how people are different) rather than general trends in behavior (i.e. tring to characterize how various factors affect decision-making within a person).
excluded = which(is.infinite(subjectData$modelAIC) | is.infinite(altSubjectData$modelAlt1AIC) | is.infinite(altSubjectData$modelAlt2AIC))
averageAIC = c(mean(subjectData$modelAIC[-excluded]), mean(altSubjectData$modelAlt1AIC[-excluded]), mean(altSubjectData$modelAlt2AIC[-excluded]))
fullAIC = length(trialData$SubjectID) * log(sum(subjectData$modelSS)/length(trialData$SubjectID)) + (2 * k * length(subjectData$SubjectID))
fullAICAlt1 = length(trialData$SubjectID) * log(sum(altSubjectData$modelAlt1SS)/length(trialData$SubjectID)) + (2 * 0 * length(subjectData$SubjectID))
fullAICAlt2 = length(trialData$SubjectID) * log(sum(altSubjectData$modelAlt2SS)/length(trialData$SubjectID)) + (2 * 0 * length(subjectData$SubjectID))
bestModel = c("Favored Model", "Alternative Model 1", "Alternative Model 2")[which(averageAIC == min(averageAIC))] #best model based on average performance per subject
bestModelFullDataset = c("Favored Model", "Alternative Model 1", "Alternative Model 2")[which(c(fullAIC, fullAICAlt1, fullAICAlt2) == min(c(fullAIC, fullAICAlt1, fullAICAlt2)))] #best model based on all data observations
excluded = find(isinf(subjectData.modelAIC) | isinf(altSubjectData.modelAlt1AIC) | isinf(altSubjectData.modelAlt2AIC));
averageAIC = [mean(subjectData.modelAIC(~excluded)), mean(altSubjectData.modelAlt1AIC(~excluded)), mean(altSubjectData.modelAlt2AIC(~excluded))];
fullAIC = length(trialData.SubjectID) * log(sum(subjectData.modelSS)/length(trialData.SubjectID)) + (2 * k * length(subjectData.SubjectID));
fullAICAlt1 = length(trialData.SubjectID) * log(sum(altSubjectData.modelAlt1SS)/length(trialData.SubjectID)) + (2 * 0 * length(subjectData.SubjectID));
fullAICAlt2 = length(trialData.SubjectID) * log(sum(altSubjectData.modelAlt2SS)/length(trialData.SubjectID)) + (2 * 0 * length(subjectData.SubjectID));
[~, bestModel] = min(averageAIC); % best model based on average performance per subject
[~, bestModelFullDataset] = min([fullAIC, fullAICAlt1, fullAICAlt2]); % best model based on all data observations
bestModel = {'Favored Model', 'Alternative Model 1', 'Alternative Model 2'}{bestModel};
bestModelFullDataset = {'Favored Model', 'Alternative Model 1', 'Alternative Model 2'}{bestModelFullDataset};
excluded = np.where(np.isinf(subjectData['modelAIC']) | np.isinf(altSubjectData['modelAlt1AIC']) | np.isinf(altSubjectData['modelAlt2AIC']))[0]
averageAIC = [np.mean(subjectData['modelAIC'][~excluded]), np.mean(altSubjectData['modelAlt1AIC'][~excluded]), np.mean(altSubjectData['modelAlt2AIC'][~excluded])]
fullAIC = len(trialData['SubjectID']) * np.log(np.sum(subjectData['modelSS']) / len(trialData['SubjectID'])) + (2 * k * len(subjectData['SubjectID']))
fullAICAlt1 = len(trialData['SubjectID']) * np.log(np.sum(altSubjectData['modelAlt1SS']) / len(trialData['SubjectID'])) + (2 * 0 * len(subjectData['SubjectID']))
fullAICAlt2 = len(trialData['SubjectID']) * np.log(np.sum(altSubjectData['modelAlt2SS']) / len(trialData['SubjectID'])) + (2 * 0 * len(subjectData['SubjectID']))
bestModel = ['Favored Model', 'Alternative Model 1', 'Alternative Model 2'][np.argmin(averageAIC)] # best model based on average performance per subject
bestModelFullDataset = ['Favored Model', 'Alternative Model 1', 'Alternative Model 2'][np.argmin([fullAIC, fullAICAlt1, fullAICAlt2])] # best model based on all data observations
Tutorials
Tutorial 1 - van Baar, Chang, & Sanfey, 2019
Create Utility Equations for our Alternative Models
utility_greed = function(greed){
return(greed)
}
utility_guilt = function(theta, greed, guilt){
return(theta * greed - (1 - theta) * guilt)
}
utility_inequity = function(theta, greed, inequity){
return(theta * greed - (1 - theta) * inequity)
}
function value = utility_greed(greed)
value = greed;
end
function value = utility_guilt(theta, greed, guilt)
value = theta * greed - (1 - theta) * guilt;
end
function value = utility_inequity(theta, greed, inequity)
value = theta * greed - (1 - theta) * inequity;
end
def utility_greed(greed):
return greed
def utility_guilt(theta, greed, guilt):
return theta * greed - (1 - theta) * guilt
def utility_inequity(theta, greed, inequity):
return theta * greed - (1 - theta) * inequity
Preallocating and Defining Functions
obj_function_guilt = function(params, df, method = "OLS") {
Theta = params[1]
predicted_utility = vector('numeric', length(df[,1]))
observed_utility = vector('numeric', length(df[,1]))
chosen = as.numeric(df[,4]) + 1
for (k in 1:length(df[,1])){
I = df[k, 2]
M = df[k, 3]
B = 4
E = 10
if (I > 10) {Choices = seq(0, (I*M), round((I*M)/10))} else {Choices = seq(0, (I * M), 1)}
Utility = vector('numeric', length(Choices))
for (n in 1:length(Choices)){
Utility[n] = utility_guilt(theta = Theta,
greed = payout_maximization(I, M, Choices[n]),
guilt = guilt(I, B, Choices[n], M))
}
predicted_utility[k] = max(Utility)
observed_utility[k] = Utility[chosen[k]]
}
if (method == "OLS"){
return(sum((predicted_utility - observed_utility)**2))
} else if (method == "MLE"){
return(-1 * sum(dnorm(observed_utility, mean = predicted_utility, sd = sd, log = TRUE)))
}
}
obj_function_inequity = function(params, df, method = "OLS") {
Theta = params[1]
predicted_utility = vector('numeric', length(df[,1]))
observed_utility = vector('numeric', length(df[,1]))
chosen = as.numeric(df[,4]) + 1
for (k in 1:length(df[,1])){
I = df[k, 2]
M = df[k, 3]
B = 4
E = 10
if (I > 10) {Choices = seq(0, (I*M), round((I*M)/10))} else {Choices = seq(0, (I * M), 1)}
Utility = vector('numeric', length(Choices))
for (n in 1:length(Choices)){
Utility[n] = utility_inequity(theta = Theta,
greed = payout_maximization(I, M, Choices[n]),
inequity = inequity(I, M, Choices[n], E))
}
predicted_utility[k] = max(Utility)
observed_utility[k] = Utility[chosen[k]]
}
if (method == "OLS"){
return(sum((predicted_utility - observed_utility)**2))
} else if (method == "MLE"){
return(-1 * sum(dnorm(observed_utility, mean = predicted_utility, sd = sd, log = TRUE)))
}
}
altSubjectData = data.frame()
function result = obj_function_guilt(params, df, method)
Theta = params(1);
predicted_utility = zeros(length(df), 1);
observed_utility = zeros(length(df), 1);
chosen = df(:, 4) + 1;
for k = 1:length(df)
I = df(k, 2);
M = df(k, 3);
B = 4;
E = 10;
if I > 10
Choices = 0:round((I * M) / 10):(I * M);
else
Choices = 0:1:(I * M);
end
Utility = zeros(size(Choices));
for n = 1:length(Choices)
Utility(n) = utility_guilt(Theta, payout_maximization(I, M, Choices(n)), guilt(I, B, Choices(n), M));
end
predicted_utility(k) = max(Utility);
observed_utility(k) = Utility(chosen(k));
end
if strcmp(method, 'OLS')
result = sum((predicted_utility - observed_utility).^2);
elseif strcmp(method, 'MLE')
result = -1 * sum(log(normpdf(observed_utility, predicted_utility, sd)));
end
end
function result = obj_function_inequity(params, df, method)
Theta = params(1);
predicted_utility = zeros(length(df), 1);
observed_utility = zeros(length(df), 1);
chosen = df(:, 4) + 1;
for k = 1:length(df)
I = df(k, 2);
M = df(k, 3);
B = 4;
E = 10;
if I > 10
Choices = 0:round((I * M) / 10):(I * M);
else
Choices = 0:1:(I * M);
end
Utility = zeros(size(Choices));
for n = 1:length(Choices)
Utility(n) = utility_inequity(Theta, payout_maximization(I, M, Choices(n)), inequity(I, M, Choices(n), E));
end
predicted_utility(k) = max(Utility);
observed_utility(k) = Utility(chosen(k));
end
if strcmp(method, 'OLS')
result = sum((predicted_utility - observed_utility).^2);
elseif strcmp(method, 'MLE')
result = -1 * sum(log(normpdf(observed_utility, predicted_utility, sd)));
end
end
altSubjectData = table();
def obj_function_guilt(params, df, method='OLS'):
Theta = params[0]
predicted_utility = np.zeros(len(df))
observed_utility = np.zeros(len(df))
chosen = df[:, 3].astype(int)
for k in range(len(df)):
I = df[k, 1]
M = df[k, 2]
B = 4
E = 10
if I > 10:
Choices = np.arange(0, (I * M) + 1, round((I * M) / 10))
else:
Choices = np.arange(0, (I * M) + 1, 1)
Utility = np.zeros(len(Choices))
for n in range(len(Choices)):
Utility[n] = utility_guilt(Theta, payout_maximization(I, M, Choices[n]), guilt(I, B, Choices[n], M))
predicted_utility[k] = np.max(Utility)
observed_utility[k] = Utility[chosen[k] - 1]
if method == 'OLS':
result = np.sum((predicted_utility - observed_utility)**2)
elif method == 'MLE':
result = -1 * np.sum(np.log(norm.pdf(observed_utility, loc=predicted_utility, scale=sd)))
return result
def obj_function_inequity(params, df, method='OLS'):
Theta = params[0]
predicted_utility = np.zeros(len(df))
observed_utility = np.zeros(len(df))
chosen = df[:, 3].astype(int)
for k in range(len(df)):
I = df[k, 1]
M = df[k, 2]
B = 4
E = 10
if I > 10:
Choices = np.arange(0, (I * M) + 1, round((I * M) / 10))
else:
Choices = np.arange(0, (I * M) + 1, 1)
Utility = np.zeros(len(Choices))
for n in range(len(Choices)):
Utility[n] = utility_inequity(Theta, payout_maximization(I, M, Choices[n]), inequity(I, M, Choices[n], E))
predicted_utility[k] = np.max(Utility)
observed_utility[k] = Utility[chosen[k] - 1]
if method == 'OLS':
result = np.sum((predicted_utility - observed_utility)**2)
elif method == 'MLE':
result = -1 * np.sum(np.log(norm.pdf(observed_utility, loc=predicted_utility, scale=sd)))
return result
altSubjectData = pd.DataFrame()
Recover Free Parameters for Each Alternative Model, Per Subject
for (i in 1:length(included_subjects)){
df = trialData[which(included_subjects[i] == trialData$Subject), ]
result_guilt = optim(fn = obj_function_guilt, par = 0.5, lower = 0, upper = 1, df = df, method = 'L-BFGS-B')
result_inequity = optim(fn = obj_function_inequity, par = 0.5, lower = 0, upper = 1, df = df, method = 'L-BFGS-B')
#Determine Model Predictions using result_guilt$pars and result_inequity$pars
}
for i = 1:length(included_subjects)
subject = included_subjects(i);
df = trialData(trialData.Subject == subject, :);
% Optimize for guilt
result_guilt = fmincon(@(theta) obj_function_guilt(theta, df), 0.5, [], [], [], [], 0, 1);
% Optimize for inequity
result_inequity = fmincon(@(theta) obj_function_inequity(theta, df), 0.5, [], [], [], [], 0, 1);
% Determine Model Predictions using result_guilt.pars and result_inequity.pars
end
for i = 1:length(included_subjects)
subject = included_subjects(i);
df = trialData(trialData.Subject == subject, :);
% Optimize for guilt
result_guilt = fmincon(@(theta) obj_function_guilt(theta, df), 0.5, [], [], [], [], 0, 1);
% Optimize for inequity
result_inequity = fmincon(@(theta) obj_function_inequity(theta, df), 0.5, [], [], [], [], 0, 1);
end
Determine Predicted Decisions for Each Alternative Model, Per Subject
for (i in 1:length(included_subjects)){
df = trialData[which(included_subjects[i] == trialData$Subject), ]
result_guilt = optim(fn = obj_function_guilt, par = 0.5, lower = 0, upper = 1, df = df, method = 'L-BFGS-B')
result_inequity = optim(fn = obj_function_inequity, par = 0.5, lower = 0, upper = 1, df = df, method = 'L-BFGS-B')
#Just Added
df$PredictionGreed = vector('numeric', length(df$Subject))
df$PredictionGuilt = df$PredictionAlt1
df$PredictionInequity = df$PredictionAlt1
for (k in 1:length(df$Returned)){
I = df$Investment[k]
M = df$Multiplier[k]
R = df$Returned[k]
B = 4
E = 10
if (I > 10) {Choices = seq(0, (I*M), round((I*M)/10))} else {Choices = seq(0, (I * M), 1)}
UtilityGreed = vector('numeric', length(Choices))
UtilityGuilt = vector('numeric', length(Choices))
UtilityInequity = vector('numeric', length(Choices))
for (n in 1:length(Choices)){
UtilityGreed[n] = utility_greed(greed = payout_maximization(I, M, Choices[n]))
UtilityGuilt[n] = utility_guilt(theta = result_guilt$par[1],
greed = payout_maximization(I, M, Choices[n]),
guilt = guilt(I, B, Choices[n], M))
UtilityInequity[n] = utility_inequity(theta = result_guilt$par[1],
greed = payout_maximization(I, M, Choices[n]),
inequity = inequity(I, M, Choices[n], E))
}
correct_choice_greed = which(UtilityGreed == max(UtilityGreed))
correct_choice_guilt = which(UtilityGuilt == max(UtilityGuilt))
correct_choice_inequity = which(UtilityInequity == max(UtilityInequity))
if (length(correct_choice_greed) > 1){
correct_choice_greed = correct_choice_greed[sample(correct_choice_greed, 1)]
}
if (length(correct_choice_guilt) > 1){
correct_choice_guilt = correct_choice_guilt[sample(correct_choice_guilt, 1)]
}
if (length(correct_choice_inequity) > 1){
correct_choice_inequity = correct_choice_inequity[sample(correct_choice_inequity, 1)]
}
df$PredictionGreed[k] = Choices[correct_choice_greed]
df$PredictionGuilt[k] = Choices[correct_choice_guilt]
df$PredictionInequity[k] = Choices[correct_choice_inequity]
}
model_NLL_Greed = -2 * log(sum(dnorm(as.numeric(df$Returned), mean = df$PredictionGreed)))
model_SS_Greed = sum((as.numeric(df$Returned) - df$PredictionGreed)**2)
model_NLL_Guilt = -2 * log(sum(dnorm(as.numeric(df$Returned), mean = df$PredictionGuilt)))
model_SS_Guilt = sum((as.numeric(df$Returned) - df$PredictionGuilt)**2)
model_NLL_Inequity = -2 * log(sum(dnorm(as.numeric(df$Returned), mean = df$PredictionInequity)))
model_SS_Inequity = sum((as.numeric(df$Returned) - df$PredictionInequity)**2)
altSubjectData[i, 1:9] = c(included_subjects[i], result_guilt$par[1], result_inequity$par[1], model_NLL_Greed, model_SS_Greed, model_NLL_Guilt, model_SS_Guilt, model_NLL_Inequity, model_SS_Inequity)
}
colnames(altSubjectData) = c('SubjectID', 'guilt_theta', 'inequity_theta', 'greed_modelNLL', 'greed_modelSS', 'guilt_modelNLL', 'guilt_modelSS', 'inequity_ModelNLL', 'inequity_ModelSS')
% Loop over included_subjects
for i = 1:length(included_subjects)
subject = included_subjects(i);
df = trialData(trialData.Subject == subject, :);
% Optimize for guilt
result_guilt = fmincon(@(theta) obj_function_guilt(theta, df), 0.5, [], [], [], [], 0, 1);
% Optimize for inequity
result_inequity = fmincon(@(theta) obj_function_inequity(theta, df), 0.5, [], [], [], [], 0, 1);
% Determine Predictions
df.PredictionGreed = df.PredictionAlt1;
df.PredictionGuilt = df.PredictionAlt1;
df.PredictionInequity = df.PredictionAlt1;
for k = 1:height(df)
I = df.Investment(k);
M = df.Multiplier(k);
R = df.Returned(k);
B = 4;
E = 10;
if I > 10
Choices = 0:round((I * M) / 10, 3):(I * M);
else
Choices = 0:1:(I * M);
end
UtilityGreed = arrayfun(@(choice) utility_greed(payout_maximization(I, M, choice)), Choices);
UtilityGuilt = arrayfun(@(choice) utility_guilt(result_guilt, payout_maximization(I, M, choice), guilt(I, B, choice, M)), Choices);
UtilityInequity = arrayfun(@(choice) utility_inequity(result_inequity, payout_maximization(I, M, choice), inequity(I, M, choice, E)), Choices);
[~, correct_choice_greed] = max(UtilityGreed);
[~, correct_choice_guilt] = max(UtilityGuilt);
[~, correct_choice_inequity] = max(UtilityInequity);
df.PredictionGreed(k) = Choices(correct_choice_greed);
df.PredictionGuilt(k) = Choices(correct_choice_guilt);
df.PredictionInequity(k) = Choices(correct_choice_inequity);
end
% Calculate model statistics
model_NLL_Greed = -2 * sum(log(normpdf(df.Returned, df.PredictionGreed, 1)));
model_SS_Greed = sum((df.Returned - df.PredictionGreed).^2);
model_NLL_Guilt = -2 * sum(log(normpdf(df.Returned, df.PredictionGuilt, 1)));
model_SS_Guilt = sum((df.Returned - df.PredictionGuilt).^2);
model_NLL_Inequity = -2 * sum(log(normpdf(df.Returned, df.PredictionInequity, 1)));
model_SS_Inequity = sum((df.Returned - df.PredictionInequity).^2);
altSubjectData(i, :) = [subject, result_guilt, result_inequity, model_NLL_Greed, model_SS_Greed, model_NLL_Guilt, model_SS_Guilt, model_NLL_Inequity, model_SS_Inequity];
end
% Set column names
colnames_altSubjectData = {'SubjectID', 'guilt_theta', 'inequity_theta', 'greed_modelNLL', 'greed_modelSS', 'guilt_modelNLL', 'guilt_modelSS', 'inequity_ModelNLL', 'inequity_ModelSS'};
altSubjectData = array2table(altSubjectData, 'VariableNames', colnames_altSubjectData);
for i, subject in enumerate(included_subjects):
df = trialData[trialData['Subject'] == subject]
# Optimize for guilt
result_guilt = minimize(obj_function_guilt, x0=[0.5], args=(df,), bounds=[(0, 1)], method='L-BFGS-B')
# Optimize for inequity
result_inequity = minimize(obj_function_inequity, x0=[0.5], args=(df,), bounds=[(0, 1)], method='L-BFGS-B')
# Determine Predictions
df['PredictionGreed'] = df['PredictionAlt1']
df['PredictionGuilt'] = df['PredictionAlt1']
df['PredictionInequity'] = df['PredictionAlt1']
for k in range(len(df['Returned'])):
I, M, R = df.loc[k, ['Investment', 'Multiplier', 'Returned']]
B, E = 4, 10
Choices = np.arange(0, I * M + 0.001, round((I * M) / 10, 3)) if I > 10 else np.arange(0, I * M + 0.001, 1)
UtilityGreed = np.array([utility_greed(payout_maximization(I, M, choice)) for choice in Choices])
UtilityGuilt = np.array([utility_guilt(result_guilt.x[0], payout_maximization(I, M, choice), guilt(I, B, choice, M)) for choice in Choices])
UtilityInequity = np.array([utility_inequity(result_inequity.x[0], payout_maximization(I, M, choice), inequity(I, M, choice, E)) for choice in Choices])
correct_choice_greed = np.argmax(UtilityGreed)
correct_choice_guilt = np.argmax(UtilityGuilt)
correct_choice_inequity = np.argmax(UtilityInequity)
df.at[k, 'PredictionGreed'] = Choices[correct_choice_greed]
df.at[k, 'PredictionGuilt'] = Choices[correct_choice_guilt]
df.at[k, 'PredictionInequity'] = Choices[correct_choice_inequity]
# Calculate model statistics
model_NLL_Greed = -2 * np.log(np.sum(norm.pdf(df['Returned'], loc=df['PredictionGreed'], scale=1)))
model_SS_Greed = np.sum((df['Returned'] - df['PredictionGreed'])**2)
model_NLL_Guilt = -2 * np.log(np.sum(norm.pdf(df['Returned'], loc=df['PredictionGuilt'], scale=1)))
model_SS_Guilt = np.sum((df['Returned'] - df['PredictionGuilt'])**2)
model_NLL_Inequity = -2 * np.log(np.sum(norm.pdf(df['Returned'], loc=df['PredictionInequity'], scale=1)))
model_SS_Inequity = np.sum((df['Returned'] - df['PredictionInequity'])**2)
altSubjectData[i, :] = [subject, result_guilt.x[0], result_inequity.x[0], model_NLL_Greed, model_SS_Greed, model_NLL_Guilt, model_SS_Guilt, model_NLL_Inequity, model_SS_Inequity]
# Set column names
colnames_altSubjectData = ['SubjectID', 'guilt_theta', 'inequity_theta', 'greed_modelNLL', 'greed_modelSS', 'guilt_modelNLL', 'guilt_modelSS', 'inequity_ModelNLL', 'inequity_ModelSS']
altSubjectData = pd.DataFrame(altSubjectData, columns=colnames_altSubjectData)
Compute Model Fit Index for Each Subject, for Each Alternative Model
N = length(altSubjectData)
altSubjectData$modelAICGreed = N * log(altSubjectData$greed_modelSS/N) + 2*0
altSubjectData$modelAICGuilt = N * log(altSubjectData$guilt_modelSS/N) + 2*1
altSubjectData$modelAICInequity = N * log(altSubjectData$inequity_ModelSS/N) + 2*1
N = length(altSubjectData);
altSubjectData.modelAICGreed = N * log(altSubjectData.greed_modelSS / N) + 2 * 0;
altSubjectData.modelAICGuilt = N * log(altSubjectData.guilt_modelSS / N) + 2 * 1;
altSubjectData.modelAICInequity = N * log(altSubjectData.inequity_ModelSS / N) + 2 * 1;
N = len(altSubjectData)
altSubjectData['modelAICGreed'] = N * np.log(altSubjectData['greed_modelSS'] / N) + 2 * 0
altSubjectData['modelAICGuilt'] = N * np.log(altSubjectData['guilt_modelSS'] / N) + 2 * 1
altSubjectData['modelAICInequity'] = N * np.log(altSubjectData['inequity_ModelSS'] / N) + 2 * 1
Compare Model Performance
#subjects excluded from analyses not included in dataset
averageAIC = c(mean(subjectData$modelAIC), mean(altSubjectData$modelAICGreed), mean(altSubjectData$modelAICGuilt), mean(altSubjectData$modelAICInequity))
fullAIC = length(trialData$Subject) * log(sum(subjectData$modelSS)/length(trialData$Subject)) + (2 * k * length(subjectData$Subject))
fullAICGreed = length(trialData$Subject) * log(sum(altSubjectData$greed_modelSS)/length(trialData$Subject)) + (2 * 0 * length(subjectData$Subject))
fullAICGuilt = length(trialData$Subject) * log(sum(altSubjectData$guilt_modelSS)/length(trialData$Subject)) + (2 * 1 * length(subjectData$Subject))
fullAICInequity = length(trialData$Subject) * log(sum(altSubjectData$inequity_ModelSS)/length(trialData$Subject)) + (2 * 1 * length(subjectData$Subject))
bestModel = c("Moral Strategies Model", "Greed Model", "Guilt Model", "Inequity Model")[which(averageAIC == min(averageAIC))] #best model based on average performance per subject
bestModelFullDataset = c("Moral Strategies Model", "Greed Model", "Guilt Model", "Inequity Model")[which(c(fullAIC, fullAICGreed, fullAICGuilt, fullAICInequity) == min(c(fullAIC, fullAICGreed, fullAICGuilt, fullAICInequity)))] #best model based on all data observations
averageAIC = [mean(subjectData.modelAIC), mean(altSubjectData.modelAICGreed), mean(altSubjectData.modelAICGuilt), mean(altSubjectData.modelAICInequity)];
fullAIC = length(trialData.Subject) * log(sum(subjectData.modelSS)/length(trialData.Subject)) + (2 * k * length(subjectData.Subject));
fullAICGreed = length(trialData.Subject) * log(sum(altSubjectData.greed_modelSS)/length(trialData.Subject)) + (2 * 0 * length(subjectData.Subject));
fullAICGuilt = length(trialData.Subject) * log(sum(altSubjectData.guilt_modelSS)/length(trialData.Subject)) + (2 * 1 * length(subjectData.Subject));
fullAICInequity = length(trialData.Subject) * log(sum(altSubjectData.inequity_ModelSS)/length(trialData.Subject)) + (2 * 1 * length(subjectData.Subject));
[~, bestModelIndex] = min(averageAIC);
bestModel = {'Moral Strategies Model', 'Greed Model', 'Guilt Model', 'Inequity Model'}{bestModelIndex}; % best model based on average performance per subject
[~, bestModelFullDatasetIndex] = min([fullAIC, fullAICGreed, fullAICGuilt, fullAICInequity]);
bestModelFullDataset = {'Moral Strategies Model', 'Greed Model', 'Guilt Model', 'Inequity Model'}{bestModelFullDatasetIndex}; % best model based on all data observations
averageAIC = [np.mean(subjectData['modelAIC']), np.mean(altSubjectData['modelAICGreed']), np.mean(altSubjectData['modelAICGuilt']), np.mean(altSubjectData['modelAICInequity'])]
fullAIC = len(trialData['Subject']) * np.log(np.sum(subjectData['modelSS']) / len(trialData['Subject'])) + (2 * k * len(subjectData['Subject']))
fullAICGreed = len(trialData['Subject']) * np.log(np.sum(altSubjectData['greed_modelSS']) / len(trialData['Subject'])) + (2 * 0 * len(subjectData['Subject']))
fullAICGuilt = len(trialData['Subject']) * np.log(np.sum(altSubjectData['guilt_modelSS']) / len(trialData['Subject'])) + (2 * 1 * len(subjectData['Subject']))
fullAICInequity = len(trialData['Subject']) * np.log(np.sum(altSubjectData['inequity_ModelSS']) / len(trialData['Subject'])) + (2 * 1 * len(subjectData['Subject']))
bestModelIndex = np.argmin(averageAIC)
bestModel = ['Moral Strategies Model', 'Greed Model', 'Guilt Model', 'Inequity Model'][bestModelIndex] # best model based on average performance per subject
bestModelFullDatasetIndex = np.argmin([fullAIC, fullAICGreed, fullAICGuilt, fullAICInequity])
bestModelFullDataset = ['Moral Strategies Model', 'Greed Model', 'Guilt Model', 'Inequity Model'][bestModelFullDatasetIndex] # best model based on all data observations
Tutorial 2 - Galvan & Sanfey, 2024
Create Utility Equations for our Alternative Models
utilityP = function(Payout){
return(Payout)
}
utilityL = function(Equality){
return(Equality)
}
utilityT = function(Equity){
return(Equity)
}
utilityPL = function(Payout, Equality, theta){
return(theta * Payout + (1 - theta) * Equality)
}
utilityPT = function(Payout, Equity, theta){
return(theta * Payout + (1 - theta) * Equity)
}
utilityLT = function(Equity, Equality, phi){
return(phi * Equality + (1 - phi) * Equity)
}
function u = utilityP(Payout)
u = Payout;
end
function u = utilityL(Equality)
u = Equality;
end
function u = utilityT(Equity)
u = Equity;
end
function u = utilityPL(Payout, Equality, theta)
u = theta * Payout + (1 - theta) * Equality;
end
function u = utilityPT(Payout, Equity, theta)
u = theta * Payout + (1 - theta) * Equity;
end
function u = utilityLT(Equity, Equality, phi)
u = phi * Equality + (1 - phi) * Equity;
end
def utilityP(Payout):
return Payout
def utilityL(Equality):
return Equality
def utilityT(Equity):
return Equity
def utilityPL(Payout, Equality, theta):
return theta * Payout + (1 - theta) * Equality
def utilityPT(Payout, Equity, theta):
return theta * Payout + (1 - theta) * Equity
def utilityLT(Equity, Equality, phi):
return phi * Equality + (1 - phi) * Equity
Preallocating and Defining Functions
obj_functionPL = function(params, df, method = "OLS") {
Theta = params[1]
predicted_utility = vector('numeric', length(df[,1]))
observed_utility = vector('numeric', length(df[,1]))
choices = seq(0, 1, 0.1)
for (k in 1:length(df[,1])){
Utility = vector('numeric', length(choices))
for (n in 1:length(choices)){
Utility[n] = utilityPL(theta = Theta,
Equality = equality(new_value(df[k, 1:10], choices[n]), df[k, 1:10], choices[n]),
Payout = payout(new_value(df[k, 1], choices[n]), df[k, 1], choices[n]))
}
predicted_utility[k] = max(Utility)
observed_utility[k] = Utility[which(df[k, 11] == choices)]
}
if (method == "OLS"){
return(sum((predicted_utility - observed_utility)**2))
} else if (method == "MLE"){
return(-1 * sum(dnorm(observed_utility, mean = predicted_utility, sd = sd, log = TRUE)))
}
}
obj_functionPT = function(params, df, method = "OLS") {
Theta = params[1]
predicted_utility = vector('numeric', length(df[,1]))
observed_utility = vector('numeric', length(df[,1]))
choices = seq(0, 1, 0.1)
for (k in 1:length(df[,1])){
Utility = vector('numeric', length(choices))
for (n in 1:length(choices)){
Utility[n] = utilityPT(theta = Theta,
Equity = equity(new_value(df[k, 1:10], choices[n]), df[k, 1:10], choices[n]),
Payout = payout(new_value(df[k, 1], choices[n]), df[k, 1], choices[n]))
}
predicted_utility[k] = max(Utility)
observed_utility[k] = Utility[which(df[k, 11] == choices)]
}
if (method == "OLS"){
return(sum((predicted_utility - observed_utility)**2))
} else if (method == "MLE"){
return(-1 * sum(dnorm(observed_utility, mean = predicted_utility, sd = sd, log = TRUE)))
}
}
obj_functionLT = function(params, df, method = "OLS") {
Phi = params[1]
predicted_utility = vector('numeric', length(df[,1]))
observed_utility = vector('numeric', length(df[,1]))
choices = seq(0, 1, 0.1)
for (k in 1:length(df[,1])){
Utility = vector('numeric', length(choices))
for (n in 1:length(choices)){
Utility[n] = utilityLT(phi = Phi,
Equity = equity(new_value(df[k, 1:10], choices[n]), df[k, 1:10], choices[n]),
Equality = equality(new_value(df[k, 1:10], choices[n]), df[k, 1:10], choices[n]))
}
predicted_utility[k] = max(Utility)
observed_utility[k] = Utility[which(df[k, 11] == choices)]
}
if (method == "OLS"){
return(sum((predicted_utility - observed_utility)**2))
} else if (method == "MLE"){
return(-1 * sum(dnorm(observed_utility, mean = predicted_utility, sd = sd, log = TRUE)))
}
}
altSubjectData = data.frame()
function objValue = obj_functionPL(params, df, method)
Theta = params(1);
predicted_utility = zeros(length(df(:, 1)), 1);
observed_utility = zeros(length(df(:, 1)), 1);
choices = 0:0.1:1;
for k = 1:length(df(:, 1))
Utility = zeros(length(choices), 1);
for n = 1:length(choices)
Utility(n) = utilityPL(Theta, ...
equality(new_value(df(k, 1:10), choices(n)), df(k, 1:10), choices(n)), ...
payout(new_value(df(k, 1), choices(n)), df(k, 1), choices(n)));
end
predicted_utility(k) = max(Utility);
observed_utility(k) = Utility(df(k, 11) == choices);
end
if strcmp(method, 'OLS')
objValue = sum((predicted_utility - observed_utility).^2);
elseif strcmp(method, 'MLE')
objValue = -sum(log(normpdf(observed_utility, predicted_utility, sd)));
end
end
function objValue = obj_functionPT(params, df, method)
Theta = params(1);
predicted_utility = zeros(length(df(:, 1)), 1);
observed_utility = zeros(length(df(:, 1)), 1);
choices = 0:0.1:1;
for k = 1:length(df(:, 1))
Utility = zeros(length(choices), 1);
for n = 1:length(choices)
Utility(n) = utilityPT(Theta, ...
equity(new_value(df(k, 1:10), choices(n)), df(k, 1:10), choices(n)), ...
payout(new_value(df(k, 1), choices(n)), df(k, 1), choices(n)));
end
predicted_utility(k) = max(Utility);
observed_utility(k) = Utility(df(k, 11) == choices);
end
if strcmp(method, 'OLS')
objValue = sum((predicted_utility - observed_utility).^2);
elseif strcmp(method, 'MLE')
objValue = -sum(log(normpdf(observed_utility, predicted_utility, sd)));
end
end
function objValue = obj_functionLT(params, df, method)
Phi = params(1);
predicted_utility = zeros(length(df(:, 1)), 1);
observed_utility = zeros(length(df(:, 1)), 1);
choices = 0:0.1:1;
for k = 1:length(df(:, 1))
Utility = zeros(length(choices), 1);
for n = 1:length(choices)
Utility(n) = utilityLT(Phi, ...
equity(new_value(df(k, 1:10), choices(n)), df(k, 1:10), choices(n)), ...
equality(new_value(df(k, 1:10), choices(n)), df(k, 1:10), choices(n)));
end
predicted_utility(k) = max(Utility);
observed_utility(k) = Utility(df(k, 11) == choices);
end
if strcmp(method, 'OLS')
objValue = sum((predicted_utility - observed_utility).^2);
elseif strcmp(method, 'MLE')
objValue = -sum(log(normpdf(observed_utility, predicted_utility, sd)));
end
end
altSubjectData = table();
def obj_functionPL(params, df, method="OLS"):
Theta = params[0]
predicted_utility = np.zeros(len(df), dtype=float)
observed_utility = np.zeros(len(df), dtype=float)
choices = np.arange(0, 1.1, 0.1)
for k in range(len(df)):
Utility = np.zeros(len(choices), dtype=float)
for n in range(len(choices)):
Utility[n] = utilityPL(Theta,
equality(new_value(df.iloc[k, 0:9], choices[n]), df.iloc[k, 0:9], choices[n]),
payout(new_value(df.iloc[k, 0], choices[n]), df.iloc[k, 0], choices[n]))
predicted_utility[k] = max(Utility)
observed_utility[k] = Utility[df.iloc[k, 10] == choices]
if method == "OLS":
obj_value = np.sum((predicted_utility - observed_utility)**2)
elif method == "MLE":
obj_value = -np.sum(np.log(norm.pdf(observed_utility, loc=predicted_utility, scale=sd)))
return obj_value
def obj_functionPT(params, df, method="OLS"):
Theta = params[0]
predicted_utility = np.zeros(len(df), dtype=float)
observed_utility = np.zeros(len(df), dtype=float)
choices = np.arange(0, 1.1, 0.1)
for k in range(len(df)):
Utility = np.zeros(len(choices), dtype=float)
for n in range(len(choices)):
Utility[n] = utilityPT(Theta,
equity(new_value(df.iloc[k, 0:9], choices[n]), df.iloc[k, 0:9], choices[n]),
payout(new_value(df.iloc[k, 0], choices[n]), df.iloc[k, 0], choices[n]))
predicted_utility[k] = max(Utility)
observed_utility[k] = Utility[df.iloc[k, 10] == choices]
if method == "OLS":
obj_value = np.sum((predicted_utility - observed_utility)**2)
elif method == "MLE":
obj_value = -np.sum(np.log(norm.pdf(observed_utility, loc=predicted_utility, scale=sd)))
return obj_value
def obj_functionLT(params, df, method="OLS"):
Phi = params[0]
predicted_utility = np.zeros(len(df), dtype=float)
observed_utility = np.zeros(len(df), dtype=float)
choices = np.arange(0, 1.1, 0.1)
for k in range(len(df)):
Utility = np.zeros(len(choices), dtype=float)
for n in range(len(choices)):
Utility[n] = utilityLT(Phi,
equity(new_value(df.iloc[k, 0:9], choices[n]), df.iloc[k, 0:9], choices[n]),
equality(new_value(df.iloc[k, 0:9], choices[n]), df.iloc[k, 0:9], choices[n]))
predicted_utility[k] = max(Utility)
observed_utility[k] = Utility[df.iloc[k, 10] == choices]
if method == "OLS":
obj_value = np.sum((predicted_utility - observed_utility)**2)
elif method == "MLE":
obj_value = -np.sum(np.log(norm.pdf(observed_utility, loc=predicted_utility, scale=sd)))
return obj_value
altSubjectData = pd.DataFrame()
Recover Free Parameters for Each Alternative Model, Per Subject Per Condition
for (i in 1:length(included_subjects)){
datafile = paste(parentfolder, included_subjects[i], restoffilepath, sep = '') # produces a character vector 'parentfolder/included_subjects[i]**.filetype'
fullData = read.csv2(datafile)
parametersPerCondition = vector('numeric', length(conditions) * 3)
SSPerCondition = vector('numeric', length(conditions) * 6)
for (j in 1:length(conditions)){
df = fullData[which(fullData$condition == conditions[j]), c(49, 40:48), 33] #49 is subject's initial allocation, 40:48 are players 1:9 initial allocation, 33 is redistribution rate
df$redistributionRate = df$redistributionRate/100 #converting to a decimal from a percent
resultPL = optim(obj_functionPL, par = 0.5, lower = 0, upper = 1, df = df, method = 'L-BFGS-B')
resultPT = optim(obj_functionPT, par = 0.5, lower = 0, upper = 1, df = df, method = 'L-BFGS-B')
resultLT = optim(obj_functionLT, par = 0.5, lower = 0, upper = 1, df = df, method = 'L-BFGS-B')
parametersPerCondition[(((j - 1) * 3) + 1):(j * 3)] = c(resultPL$par[1], resultPT$par[1], resultLT$par[1])
#Determine Predictions
}
}
parametersPerCondition = zeros(length(conditions) * 3, 1);
SSPerCondition = zeros(length(conditions) * 6, 1);
for j = 1:length(conditions)
df = fullData(fullData.condition == conditions{j}, [49, 40:48, 33]); % 49 is subject's initial allocation, 40:48 are players 1:9 initial allocation, 33 is redistribution rate
df.redistributionRate = df.redistributionRate / 100; % converting to a decimal from a percent
resultPL = fmincon(@(params) obj_functionPL(params, df, 'OLS'), 0.5, [], [], [], [], 0, 1);
resultPT = fmincon(@(params) obj_functionPT(params, df, 'OLS'), 0.5, [], [], [], [], 0, 1);
resultLT = fmincon(@(params) obj_functionLT(params, df, 'OLS'), 0.5, [], [], [], [], 0, 1);
parametersPerCondition(((j - 1) * 3) + 1 : j * 3) = [resultPL; resultPT; resultLT];
% Determine Predictions
end
from scipy.optimize import minimize
parametersPerCondition = np.zeros(len(conditions) * 3)
SSPerCondition = np.zeros(len(conditions) * 6)
for j in range(len(conditions)):
df = fullData[fullData['condition'] == conditions[j]][[47, * range(39, 48), 32]] # 49 is subject's initial allocation, 40:48 are players 1:9 initial allocation, 33 is redistribution rate
df['redistributionRate'] = df['redistributionRate'] / 100 # converting to a decimal from a percent
resultPL = minimize(lambda params: obj_functionPL(params, df, 'OLS'), 0.5, bounds=[(0, 1)])
resultPT = minimize(lambda params: obj_functionPT(params, df, 'OLS'), 0.5, bounds=[(0, 1)])
resultLT = minimize(lambda params: obj_functionLT(params, df, 'OLS'), 0.5, bounds=[(0, 1)])
parametersPerCondition[((j - 1) * 3):(j * 3)] = [resultPL.x[0], resultPT.x[0], resultLT.x[0]]
# Determine Predictions
Determine Predicted Decisions for Each Alternative Model, Per Subject Per Condition
for (i in 1:length(included_subjects)){
datafile = paste(parentfolder, included_subjects[i], restoffilepath, sep = '') # produces a character vector 'parentfolder/included_subjects[i]**.filetype'
fullData = read.csv2(datafile)
parametersPerCondition = vector('numeric', length(conditions) * 3)
SSPerCondition = vector('numeric', length(conditions) * 6)
for (j in 1:length(conditions)){
df = fullData[which(fullData$condition == conditions[j]), c(49, 40:48), 33] #49 is subject's initial allocation, 40:48 are players 1:9 initial allocation, 33 is redistribution rate
df$redistributionRate = df$redistributionRate/100 #converting to a decimal from a percent
resultPL = optim(obj_functionPL, par = 0.5, lower = 0, upper = 1, df = df, method = 'L-BFGS-B')
resultPT = optim(obj_functionPT, par = 0.5, lower = 0, upper = 1, df = df, method = 'L-BFGS-B')
resultLT = optim(obj_functionLT, par = 0.5, lower = 0, upper = 1, df = df, method = 'L-BFGS-B')
parametersPerCondition[(((j - 1) * 3) + 1):(j * 3)] = c(resultPL$par[1], resultPT$par[1], resultLT$par[1])
#Just Added
df$PredictionP = vector('numeric')
df$PredictionL = vector('numeric')
df$PredictionT = vector('numeric')
df$PredictionPT = vector('numeric')
df$PredictionPL = vector('numeric')
df$PredictionLT = vector('numeric')
for (k in 1:length(df$Decisions)){
UtilityP = vector('numeric', length(Choices))
UtilityL = vector('numeric', length(Choices))
UtilityT = vector('numeric', length(Choices))
UtilityPT = vector('numeric', length(Choices))
UtilityPL = vector('numeric', length(Choices))
UtilityLT = vector('numeric', length(Choices))
for (n in 1:length(Choices)){
UtilityP[n] = utilityP(Payout = payout(new_value(df[k, 1], choices[n]), df[k, 1], choices[n]))
UtilityL[n] = utilityL(Equality = equality(new_value(df[k, 1:10], choices[n]), df[k, 1:10], choices[n]))
UtilityT[n] = utilityT(Equity = equity(new_value(df[k, 1:10], choices[n]), df[k, 1:10], choices[n]))
UtilityPT[n] = utilityPT(theta = resultPT$par[1],
Equity = equity(new_value(df[k, 1:10], choices[n]), df[k, 1:10], choices[n]),
Payout = payout(new_value(df[k, 1], choices[n]), df[k, 1], choices[n]))
UtilityPL[n] = utilityPL(theta = resultPL$par[1],
Equality = equality(new_value(df[k, 1:10], choices[n]), df[k, 1:10], choices[n]),
Payout = payout(new_value(df[k, 1], choices[n]), df[k, 1], choices[n]))
UtilityLT[n] = utilityLT(phi = resultLT$par[1],
Equity = equity(new_value(df[k, 1:10], choices[n]), df[k, 1:10], choices[n]),
Equality = equality(new_value(df[k, 1:10], choices[n]), df[k, 1:10], choices[n]))
}
correct_choiceP = which(UtilityP == max(UtilityP))
correct_choiceL = which(UtilityL == max(UtilityL))
correct_choiceT = which(UtilityT == max(UtilityT))
correct_choicePT = which(UtilityPT == max(UtilityPT))
correct_choicePL = which(UtilityPL == max(UtilityPL))
correct_choiceLT = which(UtilityLT == max(UtilityLT))
df$predictedP[k] = new_value(df$myself[k], Choices[correct_choiceP[sample(length(correct_choiceP), 1)]])
df$predictedL[k] = new_value(df$myself[k], Choices[correct_choiceL[sample(length(correct_choiceL), 1)]])
df$predictedT[k] = new_value(df$myself[k], Choices[correct_choiceT[sample(length(correct_choiceT), 1)]])
df$predictedPL[k] = new_value(df$myself[k], Choices[correct_choicePT[sample(length(correct_choicePT), 1)]])
df$predictedPT[k] = new_value(df$myself[k], Choices[correct_choicePL[sample(length(correct_choiceTL), 1)]])
df$predictedLT[k] = new_value(df$myself[k], Choices[correct_choiceLT[sample(length(correct_choiceLT), 1)]])
}
df$Outcome = new_value(df$myself, df$redistributionRate)
SSPerCondition[(((j - 1) * 6) + 1):(j*6)] = c(sum((df$Outcome - df$predictedP)**2),
sum((df$Outcome - df$predictedL)**2),
sum((df$Outcome - df$predictedT)**2),
sum((df$Outcome - df$predictedPL)**2),
sum((df$Outcome - df$predictedPT)**2),
sum((df$Outcome - df$predictedLT)**2))
}
altSubjectData[i, 1:37] = c(included_subjects[i], parametersPerCondition, SSPerCondition)
}
colnames(altSubjectData) = c('SubjectID',
'thetaPLMerit', 'thetaPTMerit', 'phiLTMerit',
'thetaPLEntitlement', 'thetaPTEntitlement', 'phiLTEntitlement',
'thetaPLCorruption', 'thetaPTCorruption', 'phiLTCorruption',
'thetaPLLuck', 'thetaPTLuck', 'phiLTLuck',
'SSPMerit', 'SSLMerit', 'SSTMerit', 'SSPLMerit', 'SSPTMerit', 'SSLTMerit',
'SSPEntitlement', 'SSLEntitlement', 'SSTEntitlement', 'SSPLEntitlement', 'SSPTEntitlement', 'SSLTEntitlement',
'SSPCorruption', 'SSLCorruption', 'SSTCorruption', 'SSPLCorruption', 'SSPTCorruption', 'SSLTCorruption',
'SSPLuck', 'SSLLuck', 'SSTLuck', 'SSPLLuck', 'SSPLuck', 'SSLTLuck')
for i = 1:length(included_subjects)
datafile = strcat(parentfolder, included_subjects{i}, restoffilepath);
fullData = readtable(datafile);
parametersPerCondition = zeros(length(conditions) * 3, 1);
SSPerCondition = zeros(length(conditions) * 6, 1);
for j = 1:length(conditions)
df = fullData(fullData.condition == conditions(j), [49, 40:48, 33]);
df.redistributionRate = df.redistributionRate / 100;
resultPL = fmincon(@(params) obj_functionPL(params, df), 0.5, [], [], [], [], 0, 1);
resultPT = fmincon(@(params) obj_functionPT(params, df), 0.5, [], [], [], [], 0, 1);
resultLT = fmincon(@(params) obj_functionLT(params, df), 0.5, [], [], [], [], 0, 1);
parametersPerCondition(((j - 1) * 3) + 1:(j * 3)) = [resultPL(1), resultPT(1), resultLT(1)];
df.PredictionP = zeros(height(df), 1);
df.PredictionL = zeros(height(df), 1);
df.PredictionT = zeros(height(df), 1);
df.PredictionPT = zeros(height(df), 1);
df.PredictionPL = zeros(height(df), 1);
df.PredictionLT = zeros(height(df), 1);
for k = 1:length(df.Decisions)
UtilityP = zeros(length(Choices), 1);
UtilityL = zeros(length(Choices), 1);
UtilityT = zeros(length(Choices), 1);
UtilityPT = zeros(length(Choices), 1);
UtilityPL = zeros(length(Choices), 1);
UtilityLT = zeros(length(Choices), 1);
for n = 1:length(Choices)
UtilityP(n) = utilityP(payout(new_value(df{k, 1}, Choices{n}), df{k, 1}, Choices{n}));
UtilityL(n) = utilityL(equality(new_value(df{k, 1:10}, Choices{n}), df{k, 1:10}, Choices{n}));
UtilityT(n) = utilityT(equity(new_value(df{k, 1:10}, Choices{n}), df{k, 1:10}, Choices{n}));
UtilityPT(n) = utilityPT(resultPT(1), equity(new_value(df{k, 1:10}, Choices{n}), df{k, 1:10}, Choices{n}), ...
payout(new_value(df{k, 1}, Choices{n}), df{k, 1}, Choices{n}));
UtilityPL(n) = utilityPL(resultPL(1), equality(new_value(df{k, 1:10}, Choices{n}), df{k, 1:10}, Choices{n}), ...
payout(new_value(df{k, 1}, Choices{n}), df{k, 1}, Choices{n}));
UtilityLT(n) = utilityLT(resultLT(1), equity(new_value(df{k, 1:10}, Choices{n}), df{k, 1:10}, Choices{n}), ...
equality(new_value(df{k, 1:10}, Choices{n}), df{k, 1:10}, Choices{n}));
end
[~, correct_choiceP] = max(UtilityP);
[~, correct_choiceL] = max(UtilityL);
[~, correct_choiceT] = max(UtilityT);
[~, correct_choicePT] = max(UtilityPT);
[~, correct_choicePL] = max(UtilityPL);
[~, correct_choiceLT] = max(UtilityLT);
df.PredictionP(k) = new_value(df.myself{k}, Choices{correct_choiceP(randi(length(correct_choiceP)))});
df.PredictionL(k) = new_value(df.myself{k}, Choices{correct_choiceL(randi(length(correct_choiceL)))});
df.PredictionT(k) = new_value(df.myself{k}, Choices{correct_choiceT(randi(length(correct_choiceT)))});
df.PredictionPL(k) = new_value(df.myself{k}, Choices{correct_choicePT(randi(length(correct_choicePT)))});
df.PredictionPT(k) = new_value(df.myself{k}, Choices{correct_choicePL(randi(length(correct_choicePL)))});
df.PredictionLT(k) = new_value(df.myself{k}, Choices{correct_choiceLT(randi(length(correct_choiceLT)))});
end
df.Outcome = new_value(df.myself, df.redistributionRate);
SSPerCondition(((j - 1) * 6) + 1:(j * 6)) = [sum((df.Outcome - df.PredictionP).^2), ...
sum((df.Outcome - df.PredictionL).^2), ...
sum((df.Outcome - df.PredictionT).^2), ...
sum((df.Outcome - df.PredictionPL).^2), ...
sum((df.Outcome - df.PredictionPT).^2), ...
sum((df.Outcome - df.PredictionLT).^2)];
end
altSubjectData(i, 1:37) = [included_subjects{i}, parametersPerCondition', SSPerCondition'];
end
colnames_altSubjectData = {'SubjectID', 'thetaPLMerit', 'thetaPTMerit', 'phiLTMerit', 'thetaPLEntitlement', 'thetaPTEntitlement', 'phiLTEntitlement', ...
'thetaPLCorruption', 'thetaPTCorruption', 'phiLTCorruption', 'thetaPLLuck', 'thetaPTLuck', 'phiLTLuck', ...
'SSPMerit', 'SSLMerit', 'SSTMerit', 'SSPLMerit', 'SSPTMerit', 'SSLTMerit', ...
'SSPEntitlement', 'SSLEntitlement', 'SSTEntitlement', 'SSPLEntitlement', 'SSPTEntitlement', 'SSLTEntitlement', ...
'SSPCorruption', 'SSLCorruption', 'SSTCorruption', 'SSPLCorruption', 'SSPTCorruption', 'SSLTCorruption', ...
'SSPLuck', 'SSLLuck', 'SSTLuck', 'SSPLLuck', 'SSPLuck', 'SSLTLuck'};
colnames(altSubjectData) = colnames_altSubjectData;
for i in range(len(included_subjects)):
datafile = parentfolder + included_subjects[i] + restoffilepath
fullData = pd.read_csv(datafile)
parametersPerCondition = np.zeros(len(conditions) * 3)
SSPerCondition = np.zeros(len(conditions) * 6)
for j in range(len(conditions)):
df = fullData[fullData['condition'] == conditions[j]].iloc[:, [48] + list(range(39, 48)) + [32]]
df['redistributionRate'] = df['redistributionRate'] / 100
resultPL = minimize(lambda params: obj_functionPL(params, df), 0.5, bounds=[(0, 1)])
resultPT = minimize(lambda params: obj_functionPT(params, df), 0.5, bounds=[(0, 1)])
resultLT = minimize(lambda params: obj_functionLT(params, df), 0.5, bounds=[(0, 1)])
parametersPerCondition[((j - 1) * 3):(j * 3)] = [resultPL.x[0], resultPT.x[0], resultLT.x[0]]
df['PredictionP'] = np.nan
df['PredictionL'] = np.nan
df['PredictionT'] = np.nan
df['PredictionPT'] = np.nan
df['PredictionPL'] = np.nan
df['PredictionLT'] = np.nan
for k in range(len(df['Decisions'])):
UtilityP = np.zeros(len(Choices))
UtilityL = np.zeros(len(Choices))
UtilityT = np.zeros(len(Choices))
UtilityPT = np.zeros(len(Choices))
UtilityPL = np.zeros(len(Choices))
UtilityLT = np.zeros(len(Choices))
for n in range(len(Choices)):
UtilityP[n] = utilityP(Payout=payout(new_value(df.iloc[k, 0], Choices[n]), df.iloc[k, 0], Choices[n]))
UtilityL[n] = utilityL(Equality=equality(new_value(df.iloc[k, 0:9], Choices[n]), df.iloc[k, 0:9], Choices[n]))
UtilityT[n] = utilityT(Equity=equity(new_value(df.iloc[k, 0:9], Choices[n]), df.iloc[k, 0:9], Choices[n]))
UtilityPT[n] = utilityPT(theta=resultPT.x[0], Equity=equity(new_value(df.iloc[k, 0:9], Choices[n]), df.iloc[k, 0:9], Choices[n]),
Payout=payout(new_value(df.iloc[k, 0], Choices[n]), df.iloc[k, 0], Choices[n]))
UtilityPL[n] = utilityPL(theta=resultPL.x[0], Equality=equality(new_value(df.iloc[k, 0:9], Choices[n]), df.iloc[k, 0:9], Choices[n]),
Payout=payout(new_value(df.iloc[k, 0], Choices[n]), df.iloc[k, 0], Choices[n]))
UtilityLT[n] = utilityLT(phi=resultLT.x[0], Equality=equality(new_value(df.iloc[k, 0:9], Choices[n]), df.iloc[k, 0:9], Choices[n]),
Equity=equity(new_value(df.iloc[k, 0:9], Choices[n]), df.iloc[k, 0:9], Choices[n]))
correct_choiceP = np.argmax(UtilityP)
correct_choiceL = np.argmax(UtilityL)
correct_choiceT = np.argmax(UtilityT)
correct_choicePT = np.argmax(UtilityPT)
correct_choicePL = np.argmax(UtilityPL)
correct_choiceLT = np.argmax(UtilityLT)
df.at[k, 'PredictionP'] = new_value(df.iloc[k, 'myself'], Choices[correct_choiceP])
df.at[k, 'PredictionL'] = new_value(df.iloc[k, 'myself'], Choices[correct_choiceL])
df.at[k, 'PredictionT'] = new_value(df.iloc[k, 'myself'], Choices[correct_choiceT])
df.at[k, 'PredictionPL'] = new_value(df.iloc[k, 'myself'], Choices[correct_choicePT])
df.at[k, 'PredictionPT'] = new_value(df.iloc[k, 'myself'], Choices[correct_choicePL])
df.at[k, 'PredictionLT'] = new_value(df.iloc[k, 'myself'], Choices[correct_choiceLT])
df['Outcome'] = new_value(df['myself'], df['redistributionRate'])
SSPerCondition[((j - 1) * 6):(j * 6)] = [np.sum((df['Outcome'] - df['PredictionP'])**2),
np.sum((df['Outcome'] - df['PredictionL'])**2),
np.sum((df['Outcome'] - df['PredictionT'])**2),
np.sum((df['Outcome'] - df['PredictionPL'])**2),
np.sum((df['Outcome'] - df['PredictionPT'])**2),
np.sum((df['Outcome'] - df['PredictionLT'])**2)]
altSubjectData[i, 0:36] = np.concatenate(([included_subjects[i]], parametersPerCondition, SSPerCondition))
colnames = ['SubjectID', 'thetaPLMerit', 'thetaPTMerit', 'phiLTMerit', 'thetaPLEntitlement', 'thetaPTEntitlement', 'phiLTEntitlement',
'thetaPLCorruption', 'thetaPTCorruption', 'phiLTCorruption', 'thetaPLLuck', 'thetaPTLuck', 'phiLTLuck',
'SSPMerit', 'SSLMerit', 'SSTMerit', 'SSPLMerit', 'SSPTMerit', 'SSLTMerit',
'SSPEntitlement', 'SSLEntitlement', 'SSTEntitlement', 'SSPLEntitlement', 'SSPTEntitlement', 'SSLTEntitlement',
'SSPCorruption', 'SSLCorruption', 'SSTCorruption', 'SSPLCorruption', 'SSPTCorruption', 'SSLTCorruption',
'SSPLuck', 'SSLLuck', 'SSTLuck', 'SSPLLuck', 'SSPLuck', 'SSLTLuck']
altSubjectData = pd.DataFrame(altSubjectData, columns=colnames)
Compute Model Fit Index for Each Subject, for Each Alternative Model Per Condition
N = length(df[, 1])
k = 2
altSubjectData$AICPMerit = N * log(altSubjectData$SSPMerit/N) + 2*k
altSubjectData$AICLMerit = N * log(altSubjectData$SSLMerit/N) + 2*k
altSubjectData$AICTMerit = N * log(altSubjectData$SSTMerit/N) + 2*k
altSubjectData$AICPLMerit = N * log(altSubjectData$SSPLMerit/N) + 2*k
altSubjectData$AICPTMerit = N * log(altSubjectData$SSPTMerit/N) + 2*k
altSubjectData$AICLTMerit = N * log(altSubjectData$SSLTMerit/N) + 2*k
altSubjectData$AICPEntitlement = N * log(altSubjectData$SSPEntitlement/N) + 2*k
altSubjectData$AICLEntitlement = N * log(altSubjectData$SSLEntitlement/N) + 2*k
altSubjectData$AICTEntitlement = N * log(altSubjectData$SSTEntitlement/N) + 2*k
altSubjectData$AICPLEntitlement = N * log(altSubjectData$SSPLEntitlement/N) + 2*k
altSubjectData$AICPTEntitlement = N * log(altSubjectData$SSPTEntitlement/N) + 2*k
altSubjectData$AICLTEntitlement = N * log(altSubjectData$SSLTEntitlement/N) + 2*k
altSubjectData$AICPCorruption = N * log(altSubjectData$SSPCorruption/N) + 2*k
altSubjectData$AICLCorruption = N * log(altSubjectData$SSLCorruption/N) + 2*k
altSubjectData$AICTCorruption = N * log(altSubjectData$SSTCorruption/N) + 2*k
altSubjectData$AICPLCorruption = N * log(altSubjectData$SSPLCorruption/N) + 2*k
altSubjectData$AICPTCorruption = N * log(altSubjectData$SSPTCorruption/N) + 2*k
altSubjectData$AICLTCorruption = N * log(altSubjectData$SSLTCorruption/N) + 2*k
altSubjectData$AICPLuck = N * log(altSubjectData$SSPLuck/N) + 2*k
altSubjectData$AICLLuck = N * log(altSubjectData$SSLLuck/N) + 2*k
altSubjectData$AICTLuck = N * log(altSubjectData$SSTLuck/N) + 2*k
altSubjectData$AICPLLuck = N * log(altSubjectData$SSPLLuck/N) + 2*k
altSubjectData$AICPTLuck = N * log(altSubjectData$SSPTLuck/N) + 2*k
altSubjectData$AICLTLuck = N * log(altSubjectData$SSLTLuck/N) + 2*k
N = length(df(:, 1));
k = 2;
altSubjectData.AICPMerit = N * log(altSubjectData.SSPMerit / N) + 2 * k;
altSubjectData.AICLMerit = N * log(altSubjectData.SSLMerit / N) + 2 * k;
altSubjectData.AICTMerit = N * log(altSubjectData.SSTMerit / N) + 2 * k;
altSubjectData.AICPLMerit = N * log(altSubjectData.SSPLMerit / N) + 2 * k;
altSubjectData.AICPTMerit = N * log(altSubjectData.SSPTMerit / N) + 2 * k;
altSubjectData.AICLTMerit = N * log(altSubjectData.SSLTMerit / N) + 2 * k;
altSubjectData.AICPEntitlement = N * log(altSubjectData.SSPEntitlement / N) + 2 * k;
altSubjectData.AICLEntitlement = N * log(altSubjectData.SSLEntitlement / N) + 2 * k;
altSubjectData.AICTEntitlement = N * log(altSubjectData.SSTEntitlement / N) + 2 * k;
altSubjectData.AICPLEntitlement = N * log(altSubjectData.SSPLEntitlement / N) + 2 * k;
altSubjectData.AICPTEntitlement = N * log(altSubjectData.SSPTEntitlement / N) + 2 * k;
altSubjectData.AICLTEntitlement = N * log(altSubjectData.SSLTEntitlement / N) + 2 * k;
altSubjectData.AICPCorruption = N * log(altSubjectData.SSPCorruption / N) + 2 * k;
altSubjectData.AICLCorruption = N * log(altSubjectData.SSLCorruption / N) + 2 * k;
altSubjectData.AICTCorruption = N * log(altSubjectData.SSTCorruption / N) + 2 * k;
altSubjectData.AICPLCorruption = N * log(altSubjectData.SSPLCorruption / N) + 2 * k;
altSubjectData.AICPTCorruption = N * log(altSubjectData.SSPTCorruption / N) + 2 * k;
altSubjectData.AICLTCorruption = N * log(altSubjectData.SSLTCorruption / N) + 2 * k;
altSubjectData.AICPLuck = N * log(altSubjectData.SSPLuck / N) + 2 * k;
altSubjectData.AICLLuck = N * log(altSubjectData.SSLLuck / N) + 2 * k;
altSubjectData.AICTLuck = N * log(altSubjectData.SSTLuck / N) + 2 * k;
altSubjectData.AICPLLuck = N * log(altSubjectData.SSPLLuck / N) + 2 * k;
altSubjectData.AICPTLuck = N * log(altSubjectData.SSPTLuck / N) + 2 * k;
altSubjectData.AICLTLuck = N * log(altSubjectData.SSLTLuck / N) + 2 * k;
N = len(df)
k = 2
altSubjectData['AICPMerit'] = N * np.log(altSubjectData['SSPMerit'] / N) + 2 * k
altSubjectData['AICLMerit'] = N * np.log(altSubjectData['SSLMerit'] / N) + 2 * k
altSubjectData['AICTMerit'] = N * np.log(altSubjectData['SSTMerit'] / N) + 2 * k
altSubjectData['AICPLMerit'] = N * np.log(altSubjectData['SSPLMerit'] / N) + 2 * k
altSubjectData['AICPTMerit'] = N * np.log(altSubjectData['SSPTMerit'] / N) + 2 * k
altSubjectData['AICLTMerit'] = N * np.log(altSubjectData['SSLTMerit'] / N) + 2 * k
altSubjectData['AICPEntitlement'] = N * np.log(altSubjectData['SSPEntitlement'] / N) + 2 * k
altSubjectData['AICLEntitlement'] = N * np.log(altSubjectData['SSLEntitlement'] / N) + 2 * k
altSubjectData['AICTEntitlement'] = N * np.log(altSubjectData['SSTEntitlement'] / N) + 2 * k
altSubjectData['AICPLEntitlement'] = N * np.log(altSubjectData['SSPLEntitlement'] / N) + 2 * k
altSubjectData['AICPTEntitlement'] = N * np.log(altSubjectData['SSPTEntitlement'] / N) + 2 * k
altSubjectData['AICLTEntitlement'] = N * np.log(altSubjectData['SSLTEntitlement'] / N) + 2 * k
altSubjectData['AICPCorruption'] = N * np.log(altSubjectData['SSPCorruption'] / N) + 2 * k
altSubjectData['AICLCorruption'] = N * np.log(altSubjectData['SSLCorruption'] / N) + 2 * k
altSubjectData['AICTCorruption'] = N * np.log(altSubjectData['SSTCorruption'] / N) + 2 * k
altSubjectData['AICPLCorruption'] = N * np.log(altSubjectData['SSPLCorruption'] / N) + 2 * k
altSubjectData['AICPTCorruption'] = N * np.log(altSubjectData['SSPTCorruption'] / N) + 2 * k
altSubjectData['AICLTCorruption'] = N * np.log(altSubjectData['SSLTCorruption'] / N) + 2 * k
altSubjectData['AICPLuck'] = N * np.log(altSubjectData['SSPLuck'] / N) + 2 * k
altSubjectData['AICLLuck'] = N * np.log(altSubjectData['SSLLuck'] / N) + 2 * k
altSubjectData['AICTLuck'] = N * np.log(altSubjectData['SSTLuck'] / N) + 2 * k
altSubjectData['AICPLLuck'] = N * np.log(altSubjectData['SSPLLuck'] / N) + 2 * k
altSubjectData['AICPTLuck'] = N * np.log(altSubjectData['SSPTLuck'] / N) + 2 * k
altSubjectData['AICLTLuck'] = N * np.log(altSubjectData['SSLTLuck'] / N) + 2 * k
Compare Model Performance
excludedM = which(is.infinite(subjectData$AICMerit)) #any perfectly predicted subjects by alternative models should also be perfectly predicted by the most complex model as well
excludedE = which(is.infinite(subjectData$AICEntitlement))
excludedC = which(is.infinite(subjectData$AICCorruption))
excludedL = which(is.infinite(subjectData$AICLuck))
models = c('Favored Model', 'Equality Only', 'Equity Only', 'Payout Only', 'Payout-Equality', 'Payout-Equity', 'Equality-Equity')
bestModelMerit = models[which(c(mean(subjectData$AICMerit[-excluded]),
mean(altSubjectData$AICPMerit[-excluded]),
mean(altSubjectData$AICLMerit[-excluded]),
mean(altSubjectData$AICTMerit[-excluded]),
mean(altSubjectData$AICPLMerit[-excluded]),
mean(altSubjectData$AICPTMerit[-excluded]),
mean(altSubjectData$AICLTMerit[-excluded])) == max(c(mean(subjectData$AICMerit[-excluded]),
mean(altSubjectData$AICPMerit[-excluded]),
mean(altSubjectData$AICLMerit[-excluded]),
mean(altSubjectData$AICTMerit[-excluded]),
mean(altSubjectData$AICPLMerit[-excluded]),
mean(altSubjectData$AICPTMerit[-excluded]),
mean(altSubjectData$AICLTMerit[-excluded]))))]
bestModelEntitlement = models[which(c(mean(subjectData$AICEntitlement[-excluded]),
mean(altSubjectData$AICPEntitlement[-excluded]),
mean(altSubjectData$AICLEntitlement[-excluded]),
mean(altSubjectData$AICTEntitlement[-excluded]),
mean(altSubjectData$AICPLEntitlement[-excluded]),
mean(altSubjectData$AICPTEntitlement[-excluded]),
mean(altSubjectData$AICLTEntitlement[-excluded])) == max(c(mean(subjectData$AICEntitlement[-excluded]),
mean(altSubjectData$AICPEntitlement[-excluded]),
mean(altSubjectData$AICLEntitlement[-excluded]),
mean(altSubjectData$AICTEntitlement[-excluded]),
mean(altSubjectData$AICPLEntitlement[-excluded]),
mean(altSubjectData$AICPTEntitlement[-excluded]),
mean(altSubjectData$AICLTEntitlement[-excluded]))))]
bestModelCorruption = models[which(c(mean(subjectData$AICCorruption[-excluded]),
mean(altSubjectData$AICPCorruption[-excluded]),
mean(altSubjectData$AICLCorruption[-excluded]),
mean(altSubjectData$AICTCorruption[-excluded]),
mean(altSubjectData$AICPLCorruption[-excluded]),
mean(altSubjectData$AICPTCorruption[-excluded]),
mean(altSubjectData$AICLTCorruption[-excluded])) == max(c(mean(subjectData$AICCorruption[-excluded]),
mean(altSubjectData$AICPCorruption[-excluded]),
mean(altSubjectData$AICLCorruption[-excluded]),
mean(altSubjectData$AICTCorruption[-excluded]),
mean(altSubjectData$AICPLCorruption[-excluded]),
mean(altSubjectData$AICPTCorruption[-excluded]),
mean(altSubjectData$AICLTCorruption[-excluded]))))]
bestModelLuck = models[which(c(mean(subjectData$AICLuck[-excluded]),
mean(altSubjectData$AICPLuck[-excluded]),
mean(altSubjectData$AICLLuck[-excluded]),
mean(altSubjectData$AICTLuck[-excluded]),
mean(altSubjectData$AICPLLuck[-excluded]),
mean(altSubjectData$AICPTLuck[-excluded]),
mean(altSubjectData$AICLTLuck[-excluded])) == max(c(mean(subjectData$AICLuck[-excluded]),
mean(altSubjectData$AICPLuck[-excluded]),
mean(altSubjectData$AICLLuck[-excluded]),
mean(altSubjectData$AICTLuck[-excluded]),
mean(altSubjectData$AICPLLuck[-excluded]),
mean(altSubjectData$AICPTLuck[-excluded]),
mean(altSubjectData$AICLTLuck[-excluded]))))]
excludedM = find(isinf(subjectData.AICMerit)); % any perfectly predicted subjects by alternative models should also be perfectly predicted by the most complex model as well
excludedE = find(isinf(subjectData.AICEntitlement));
excludedC = find(isinf(subjectData.AICCorruption));
excludedL = find(isinf(subjectData.AICLuck));
models = {'Favored Model', 'Equality Only', 'Equity Only', 'Payout Only', 'Payout-Equality', 'Payout-Equity', 'Equality-Equity'};
meanAICMerit = mean(subjectData.AICMerit(~excludedM));
meanAICPEMerit = mean(altSubjectData.AICPMerit(~excludedM));
meanAICLEMerit = mean(altSubjectData.AICLMerit(~excludedM));
meanAICTEMerit = mean(altSubjectData.AICTMerit(~excludedM));
meanAICPLEMerit = mean(altSubjectData.AICPLMerit(~excludedM));
meanAICPTEMerit = mean(altSubjectData.AICPTMerit(~excludedM));
meanAICLTEMerit = mean(altSubjectData.AICLTMerit(~excludedM));
meanAICs = [meanAICMerit, meanAICPEMerit, meanAICLEMerit, meanAICTEMerit, meanAICPLEMerit, meanAICPTEMerit, meanAICLTEMerit];
[~, bestModelIndex] = max(meanAICs);
bestModelMerit = models{bestModelIndex};meanAICEntitlement = mean(subjectData.AICEntitlement(~excludedE));
meanAICPEEntitlement = mean(altSubjectData.AICPEntitlement(~excludedE));
meanAICLEEntitlement = mean(altSubjectData.AICLEntitlement(~excludedE));
meanAICTEEntitlement = mean(altSubjectData.AICTEntitlement(~excludedE));
meanAICPLEEntitlement = mean(altSubjectData.AICPLEntitlement(~excludedE));
meanAICPTEEntitlement = mean(altSubjectData.AICPTEntitlement(~excludedE));
meanAICLTEEntitlement = mean(altSubjectData.AICLTEntitlement(~excludedE));
meanAICsEntitlement = [meanAICEntitlement, meanAICPEEntitlement, meanAICLEEntitlement, meanAICTEEntitlement, meanAICPLEEntitlement, meanAICPTEEntitlement, meanAICLTEEntitlement];
[~, bestModelIndexEntitlement] = max(meanAICsEntitlement);
bestModelEntitlement = models{bestModelIndexEntitlement};
meanAICCorruption = mean(subjectData.AICCorruption(~excludedC));
meanAICPCorruption = mean(altSubjectData.AICPCorruption(~excludedC));
meanAICLCorruption = mean(altSubjectData.AICLCorruption(~excludedC));
meanAICTCorruption = mean(altSubjectData.AICTCorruption(~excludedC));
meanAICPLCorruption = mean(altSubjectData.AICPLCorruption(~excludedC));
meanAICPTCorruption = mean(altSubjectData.AICPTCorruption(~excludedC));
meanAICLTCorruption = mean(altSubjectData.AICLTCorruption(~excludedC));
meanAICsCorruption = [meanAICCorruption, meanAICPCorruption, meanAICLCorruption, meanAICTCorruption, meanAICPLCorruption, meanAICPTCorruption, meanAICLTCorruption];
[~, bestModelIndexCorruption] = max(meanAICsCorruption);
bestModelCorruption = models{bestModelIndexCorruption};
meanAICLuck = mean(subjectData.AICLuck(~excludedL));
meanAICPLuck = mean(altSubjectData.AICPLuck(~excludedL));
meanAICLLuck = mean(altSubjectData.AICLLuck(~excludedL));
meanAICTLuck = mean(altSubjectData.AICTLuck(~excludedL));
meanAICPLLuck = mean(altSubjectData.AICPLLuck(~excludedL));
meanAICPTLuck = mean(altSubjectData.AICPTLuck(~excludedL));
meanAICLTLuck = mean(altSubjectData.AICLTLuck(~excludedL));
meanAICsLuck = [meanAICLuck, meanAICPLuck, meanAICLLuck, meanAICTLuck, meanAICPLLuck, meanAICPTLuck, meanAICLTLuck];
[~, bestModelIndexLuck] = max(meanAICsLuck);
bestModelLuck = models{bestModelIndexLuck};
# Find indices of subjects with infinite AIC values
excludedM = np.where(np.isinf(subjectData['AICMerit']))[0]
excludedE = np.where(np.isinf(subjectData['AICEntitlement']))[0]
excludedC = np.where(np.isinf(subjectData['AICCorruption']))[0]
excludedL = np.where(np.isinf(subjectData['AICLuck']))[0]
models = ['Favored Model', 'Equality Only', 'Equity Only', 'Payout Only', 'Payout-Equality', 'Payout-Equity', 'Equality-Equity']
# Calculate means for each model and find the best model for each category
bestModelMerit = models[np.argmax([np.mean(subjectData['AICMerit'][-excludedM]),
np.mean(altSubjectData['AICPMerit'][-excludedM]),
np.mean(altSubjectData['AICLMerit'][-excludedM]),
np.mean(altSubjectData['AICTMerit'][-excludedM]),
np.mean(altSubjectData['AICPLMerit'][-excludedM]),
np.mean(altSubjectData['AICPTMerit'][-excludedM]),
np.mean(altSubjectData['AICLTMerit'][-excludedM])])]
bestModelEntitlement = models[np.argmax([np.mean(subjectData['AICEntitlement'][-excludedE]),
np.mean(altSubjectData['AICPEntitlement'][-excludedE]),
np.mean(altSubjectData['AICLEntitlement'][-excludedE]),
np.mean(altSubjectData['AICTEntitlement'][-excludedE]),
np.mean(altSubjectData['AICPLEntitlement'][-excludedE]),
np.mean(altSubjectData['AICPTEntitlement'][-excludedE]),
np.mean(altSubjectData['AICLTEntitlement'][-excludedE])])]
bestModelCorruption = models[np.argmax([np.mean(subjectData['AICCorruption'][-excludedC]),
np.mean(altSubjectData['AICPCorruption'][-excludedC]),
np.mean(altSubjectData['AICLCorruption'][-excludedC]),
np.mean(altSubjectData['AICTCorruption'][-excludedC]),
np.mean(altSubjectData['AICPLCorruption'][-excludedC]),
np.mean(altSubjectData['AICPTCorruption'][-excludedC]),
np.mean(altSubjectData['AICLTCorruption'][-excludedC])])]
bestModelLuck = models[np.argmax([np.mean(subjectData['AICLuck'][-excludedL]),
np.mean(altSubjectData['AICPLuck'][-excludedL]),
np.mean(altSubjectData['AICLLuck'][-excludedL]),
np.mean(altSubjectData['AICTLuck'][-excludedL]),
np.mean(altSubjectData['AICPLLuck'][-excludedL]),
np.mean(altSubjectData['AICPTLuck'][-excludedL]),
np.mean(altSubjectData['AICLTLuck'][-excludedL])])]
Tutorial 3 - Crockett et al., 2014
Preallocating and Defining Functions
#adjusting
optimize = function(of, using, df){
tryCatch({
fmincon(of, x0 = initial_params[using], lb = lower_bounds[using], ub = upper_bounds[using], df = df, optimMethod = "MLE")
}, error = function(e){
fmincon(of, x0 = initial_params[using], lb = lower_bounds[using], ub = upper_bounds[using], df = df, optimMethod = "OLS")
})
}
obj_functionNoGamma = function(params, df, optimMethod = "MLE") {
#we want kappa and beta per condition, so we'll assign those based on the condition
#we need to read the codebook to change the columns to reflect how the data is stored in the .mat files
Epsilon = params[5]
Gamma = 0 #0 means no bias, everything else is the same
ProbHarm = vector('numeric', length(df[,1]))
choices = as.numeric(df[, 11])
for (k in 1:length(df[,1])){
isLeft = as.logical(df[k, 8])
isSelf = as.logical(df[k, 10])
if (isSelf){Kappa = params[1]; Beta = params[3]} else {Kappa = params[2]; Beta = params[4]}
moneyMore = as.numeric(df[k, 4])
moneyLess = as.numeric(df[k, 2])
shocksMore = as.numeric(df[k, 7])
shocksLess = as.numeric(df[k, 5])
Utility = utility(Harm = harm(shocksMore, shocksLess),
Payout = payout(moneyMore, moneyLess),
kappa = Kappa)
ProbHarm[k] = max(min(probability(Beta, Epsilon, Gamma, Utility, isLeft), 0.9999999999), 0.00000000001)
}
if (optimMethod == "OLS"){
return(sum((choices - ProbHarm)**2))
} else if (optimMethod == "MLE"){
return(-sum(choices * log(ProbHarm) + (1 - choices) * log(1 - ProbHarm)))
}
}
obj_functionNoEpsilon = function(params, df, optimMethod = "MLE") {
#we want kappa and beta per condition, so we'll assign those based on the condition
#we need to read the codebook to change the columns to reflect how the data is stored in the .mat files
Epsilon = 0 #0 means no noise
Gamma = 0 #and 0 means no bias, everything else is the same
ProbHarm = vector('numeric', length(df[,1]))
choices = as.numeric(df[, 11])
for (k in 1:length(df[,1])){
isLeft = as.logical(df[k, 8])
isSelf = as.logical(df[k, 10])
if (isSelf){Kappa = params[1]; Beta = params[3]} else {Kappa = params[2]; Beta = params[4]}
moneyMore = as.numeric(df[k, 4])
moneyLess = as.numeric(df[k, 2])
shocksMore = as.numeric(df[k, 7])
shocksLess = as.numeric(df[k, 5])
Utility = utility(Harm = harm(shocksMore, shocksLess),
Payout = payout(moneyMore, moneyLess),
kappa = Kappa)
ProbHarm[k] = max(min(probability(Beta, Epsilon, Gamma, Utility, isLeft), 0.9999999999), 0.00000000001)
}
if (optimMethod == "OLS"){
return(sum((choices - ProbHarm)**2))
} else if (optimMethod == "MLE"){
return(-sum(choices * log(ProbHarm) + (1 - choices) * log(1 - ProbHarm)))
}
}
obj_functionNoConditions = function(params, df, optimMethod = "MLE") {
#we want kappa and beta per condition, so we'll assign those based on the condition
#we need to read the codebook to change the columns to reflect how the data is stored in the .mat files
Kappa = params[1] #moved up here because not changing depending on trial
Beta = params[2] #same as above
Epsilon = params[3] #change from 5 since 2 fewer parameters
Gamma = params[4] #change from 6 for same reason
ProbHarm = vector('numeric', length(df[,1]))
choices = as.numeric(df[, 11])
for (k in 1:length(df[,1])){
isLeft = as.logical(df[k, 8])
isSelf = as.logical(df[k, 10])
moneyMore = as.numeric(df[k, 4])
moneyLess = as.numeric(df[k, 2])
shocksMore = as.numeric(df[k, 7])
shocksLess = as.numeric(df[k, 5])
Utility = utility(Harm = harm(shocksMore, shocksLess),
Payout = payout(moneyMore, moneyLess),
kappa = Kappa)
ProbHarm[k] = max(min(probability(Beta, Epsilon, Gamma, Utility, isLeft), 0.9999999999), 0.00000000001)
}
if (optimMethod == "OLS"){
return(sum((choices - ProbHarm)**2))
} else if (optimMethod == "MLE"){
return(-sum(choices * log(ProbHarm) + (1 - choices) * log(1 - ProbHarm)))
}
}
obj_functionNoKappa = function(params, df, optimMethod = "MLE") {
#we want kappa and beta per condition, so we'll assign those based on the condition
#we need to read the codebook to change the columns to reflect how the data is stored in the .mat files
Kappa = 0 #nonfactor
Beta = 0 #nonfactor
Epsilon = 0.5 #the only factor is noise and bias
Gamma = params[1] #bias can vary
ProbHarm = vector('numeric', length(df[,1]))
choices = as.numeric(df[, 11])
for (k in 1:length(df[,1])){
isLeft = as.logical(df[k, 8])
isSelf = as.logical(df[k, 10])
moneyMore = as.numeric(df[k, 4])
moneyLess = as.numeric(df[k, 2])
shocksMore = as.numeric(df[k, 7])
shocksLess = as.numeric(df[k, 5])
Utility = utility(Harm = harm(shocksMore, shocksLess),
Payout = payout(moneyMore, moneyLess),
kappa = Kappa)
ProbHarm[k] = max(min(probability(Beta, Epsilon, Gamma, Utility, isLeft), 0.9999999999), 0.00000000001)
}
if (optimMethod == "OLS"){
return(sum((choices - ProbHarm)**2))
} else if (optimMethod == "MLE"){
return(-sum(choices * log(ProbHarm) + (1 - choices) * log(1 - ProbHarm)))
}
}
NoGamma = seq(1, 5)
NoEpsilon = seq(1, 4)
NoConditions = c(1, 3, 5, 6)
NoKappa = 6
altSubjectData = data.frame()
Recover Free Parameters and Determine Predicted Decisions
for (i in 1:length(included_subjects)){
df = grab_data(included_subjects[i])
if (is.character(df)){next}
resultNoGamma = optimize(obj_functionNoGamma, NoGamma, df)
resultNoEpsilon = optimize(obj_functionNoEpsilon, NoEpsilon, df)
resultNoConditions = optimize(obj_functionNoConditions, NoConditions, df)
resultNoKappa = optim(obj_functionNoKappa,par = initial_params[NoKappa], lower = lower_bounds[NoKappa], upper = upper_bounds[NoKappa], df = df, method = 'L-BFGS-B')
dfNoGamma = generate_predictions(df, data.frame(par = c(resultNoGamma$par[1:5], 0)))
dfNoEpsilon = generate_predictions(df, data.frame(par = c(resultNoEpsilon$par[1:4], 0, 0)))
dfNoConditions = generate_predictions(df, data.frame(par = c(resultNoConditions$par[c(1, 1, 2, 2, 3, 4)])))
dfNoKappa = generate_predictions(df, data.frame(par = c(0, 0, 0, 0, 0.5, resultNoKappa$par[1])))
modelNG_SS = sum((dfNoGamma$choseHarm - dfNoGamma$ProbHarm)**2)
modelNG_NLL = -2*sum(dfNoGamma$choseHarm * log(dfNoGamma$ProbHarm) + (1 - dfNoGamma$choseHarm) * log(1 - dfNoGamma$ProbHarm))
modelNE_SS = sum((dfNoEpsilon$choseHarm - dfNoEpsilon$ProbHarm)**2)
modelNE_NLL = -2*sum(dfNoEpsilon$choseHarm * log(dfNoEpsilon$ProbHarm) + (1 - dfNoEpsilon$choseHarm) * log(1 - dfNoEpsilon$ProbHarm))
modelNC_SS = sum((dfNoConditions$choseHarm - dfNoConditions$ProbHarm)**2)
modelNC_NLL = -2*sum(dfNoConditions$choseHarm * log(dfNoConditions$ProbHarm) + (1 - dfNoConditions$choseHarm) * log(1 - dfNoConditions$ProbHarm))
modelNK_SS = sum((dfNoKappa$choseHarm - dfNoKappa$ProbHarm)**2)
modelNK_NLL = -2*sum(dfNoKappa$choseHarm * log(dfNoKappa$ProbHarm) + (1 - dfNoKappa$choseHarm) * log(1 - dfNoKappa$ProbHarm))
altSubjectData[i, 1:23] = c(included_subjects[i],
resultNoGamma$par[1], resultNoGamma$par[2], resultNoGamma$par[3], resultNoGamma$par[4], resultNoGamma$par[5],
resultNoEpsilon$par[1], resultNoEpsilon$par[2], resultNoEpsilon$par[3], resultNoEpsilon$par[4],
resultNoConditions$par[1], resultNoConditions$par[2], resultNoConditions$par[3], resultNoConditions$par[4],
resultNoKappa$par[1],
modelNG_SS, modelNG_NLL, modelNE_SS, modelNE_NLL, modelNC_SS, modelNC_NLL, modelNK_SS, modelNK_NLL)
df = cbind(data.frame(SubjectID = rep(included_subjects[i], length(df$trialNumber))), df)
df$ProbHarmNG = dfNoGamma$ProbHarm
df$ProbHarmNE = dfNoEpsilon$ProbHarm
df$ProbHarmNC = dfNoConditions$ProbHarm
df$ProbHarmNK = dfNoKappa$ProbHarm
if (i == 1) {altTrialData = df} else {altTrialData = rbind(altTrialData, df)}
}
colnames(altSubjectData) = c('SubjectID', 'KappaSelfM2', 'KappaOtherM2', 'BetaSelfM2', 'BetaOtherM2', 'EpsilonM2',
'KappaSelfM3', 'KappaOtherM3', 'BetaSelfM3', 'BetaOtherM3',
'KappaM4', 'BetaM4', 'EpsilonM4', 'GammaM4',
'GammaM5',
'SSM2', 'DevianceM2', 'SSM3', 'DevianceM3', 'SSM4', 'DevianceM4', 'SSM5', 'DevianceM5')
for (i in 2:ncol(altSubjectData)){altSubjectData[,i] = as.numeric(altSubjectData[,i])}
Compute Model Fit Index for Each Subject, for Each Alternative Model
altSubjectData$BICM2 = as.numeric(altSubjectData$DevianceM2) + log(152) * 5
altSubjectData$BICM3 = as.numeric(altSubjectData$DevianceM3) + log(152) * 4
altSubjectData$BICM4 = as.numeric(altSubjectData$DevianceM4) + log(152) * 4
altSubjectData$BICM5 = as.numeric(altSubjectData$DevianceM5) + log(152) * 1
Compare Model Performance
modelBIC = c(sum(subjectData$BIC[-20]), sum(altSubjectData$BICM2[-20]), sum(altSubjectData$BICM3[-20]), sum(altSubjectData$BICM4[-20]), sum(altSubjectData$BICM5[-20]))
which(modelBIC == min(modelBIC))
Tutorial 4 - Li et al., 2022
Preallocating and Defining Functions
#adjusting
of_alphaOnly = function(params, df, optimMethod = "MLE"){
params_new = rep(0, 6); params_new[c(1, 4:6)] = params
Prob1 = generatePredictions(params_new, df); Chose1 = df[,7]
if (optimMethod == "OLS"){return(sum((Chose1 - Prob1)**2))
} else if (optimMethod == "MLE"){return(-sum(Chose1 * log(Prob1) + (1 - Chose1) * log(1 - Prob1)))}
}
of_deltaOnly = function(params, df, optimMethod = "MLE"){
params_new = rep(0, 6); params_new[c(2, 4:6)] = params
Prob1 = generatePredictions(params_new, df); Chose1 = df[,7]
if (optimMethod == "OLS"){return(sum((Chose1 - Prob1)**2))
} else if (optimMethod == "MLE"){return(-sum(Chose1 * log(Prob1) + (1 - Chose1) * log(1 - Prob1)))}
}
of_rhoOnly = function(params, df, optimMethod = "MLE"){
params_new = rep(0, 6); params_new[c(3:6)] = params
Prob1 = generatePredictions(params_new, df); Chose1 = df[,7]
if (optimMethod == "OLS"){return(sum((Chose1 - Prob1)**2))
} else if (optimMethod == "MLE"){return(-sum(Chose1 * log(Prob1) + (1 - Chose1) * log(1 - Prob1)))}
}
of_ad = function(params, df, optimMethod = "MLE"){
params_new = rep(0, 6); params_new[c(1:2, 4:6)] = params
Prob1 = generatePredictions(params_new, df); Chose1 = df[,7]
if (optimMethod == "OLS"){return(sum((Chose1 - Prob1)**2))
} else if (optimMethod == "MLE"){return(-sum(Chose1 * log(Prob1) + (1 - Chose1) * log(1 - Prob1)))}
}
of_ar = function(params, df, optimMethod = "MLE"){
params_new = rep(0, 6); params_new[c(1, 3:6)] = params
Prob1 = generatePredictions(params_new, df); Chose1 = df[,7]
if (optimMethod == "OLS"){return(sum((Chose1 - Prob1)**2))
} else if (optimMethod == "MLE"){return(-sum(Chose1 * log(Prob1) + (1 - Chose1) * log(1 - Prob1)))}
}
of_dr = function(params, df, optimMethod = "MLE"){
params_new = rep(0, 6); params_new[c(2:6)] = params
Prob1 = generatePredictions(params_new, df); Chose1 = df[,7]
if (optimMethod == "OLS"){return(sum((Chose1 - Prob1)**2))
} else if (optimMethod == "MLE"){return(-sum(Chose1 * log(Prob1) + (1 - Chose1) * log(1 - Prob1)))}
}
of_noEpsilon = function(params, df, optimMethod = "MLE"){
params_new = rep(0, 6); params_new[c(1:4)] = params
Prob1 = generatePredictions(params_new, df); Chose1 = df[,7]
if (optimMethod == "OLS"){return(sum((Chose1 - Prob1)**2))
} else if (optimMethod == "MLE"){return(-sum(Chose1 * log(Prob1) + (1 - Chose1) * log(1 - Prob1)))}
}
of_noGamma = function(params, df, optimMethod = "MLE"){
params_new = rep(0, 6); params_new[c(1:5)] = params
Prob1 = generatePredictions(params_new, df); Chose1 = df[,7]
if (optimMethod == "OLS"){return(sum((Chose1 - Prob1)**2))
} else if (optimMethod == "MLE"){return(-sum(Chose1 * log(Prob1) + (1 - Chose1) * log(1 - Prob1)))}
}
of_GammaOnly = function(params, df, optimMethod = "MLE"){
params_new = rep(0, 6); params_new[5] = 0.5; params_new[6] = params
Prob1 = generatePredictions(params_new, df); Chose1 = df[,7]
if (optimMethod == "OLS"){return(sum((Chose1 - Prob1)**2))
} else if (optimMethod == "MLE"){return(-sum(Chose1 * log(Prob1) + (1 - Chose1) * log(1 - Prob1)))}
}
ofs = list(of_alphaOnly, of_deltaOnly, of_rhoOnly, of_ad, of_ar, of_dr, of_noEpsilon, of_noGamma, of_GammaOnly)
idxs = list(c(1,4:6), c(2,4:6), c(3:6), c(1:2, 4:6), c(1, 3:6), c(2:6), c(1:4), c(1:5), c(6))
altSubjectData = data.frame()
altTrialData = trialData[,-9]
altTrialData$alphaOnly_Prob1 = 0
altTrialData$deltaOnly_Prob1 = 0
altTrialData$rhoOnly_Prob1 = 0
altTrialData$ad_Prob1 = 0
altTrialData$ar_Prob1 = 0
altTrialData$dr_Prob1 = 0
altTrialData$noEpsilon_Prob1 = 0
altTrialData$noGamma_Prob1 = 0
altTrialData$gammaOnly_Prob1 = 0
Recover Free Parameters and Determine Predicted Decisions
for (i in 1:length(included_subjects)){
df = grab_data(included_subjects[i])
outputs = vector('numeric', length = sum(length(c(c(1,4:6), c(2,4:6), c(3:6), c(1:2, 4:6), c(1, 3:6), c(2:6), c(1:4), c(1:5), c(6)))) + 2*9)
j = 0
for (k in 1:length(idxs)){
idx = idxs[k][[1]]; initials = initial_params[idx]; uppers = upper_bounds[idx]; lowers = lower_bounds[idx]; of = ofs[k][[1]]
if (length(idx) > 1){
result = optimize(of, initials, lowers, uppers, df)
} else {
result = optim(par = initials, fn = of, method = "L-BFGS-B", lower = lowers, upper = uppers, df = df)
}
pars = rep(0, times = 6)
pars[idx] = result$par
df$Prob1 = 0; df$Prob1 = generatePredictions(pars, df)
model_SS = sum((df$Chose1 - df$Prob1)**2)
model_NLL = -2*sum(df$Chose1 * log(df$Prob1) + (1 - df$Chose1) * log(1 - df$Prob1))
outputs[(j+1):(j+2+length(result$par))] = c(result$par, model_SS, model_NLL)
j = j+2+length(result$par)
altTrialData[which(altTrialData$SubjectID == included_subjects[i]), (8+k)] = df$Prob1
}
altSubjectData[i, 1:56] = c(included_subjects[i], outputs)
}
colnames(altSubjectData) = c('SubjectID',
'Alpha_M1', 'Beta_M1', 'Epsilon_M1', 'Gamma_M1', 'SS_M1', 'Deviance_M1',
'Delta_M2', 'Beta_M2', 'Epsilon_M2', 'Gamma_M2', 'SS_M2', 'Deviance_M2',
'Rho_M3', 'Beta_M3', 'Epsilon_M3', 'Gamma_M3', 'SS_M3', 'Deviance_M3',
'Alpha_M4', 'Delta_M4', 'Beta_M4', 'Epsilon_M4', 'Gamma_M4', 'SS_M4', 'Deviance_M4',
'Alpha_M5', 'Rho_M5', 'Beta_M5', 'Epsilon_M5', 'Gamma_M5', 'SS_M5', 'Deviance_M5',
'Delta_M6', 'Rho_M6', 'Beta_M6', 'Epsilon_M6', 'Gamma_M6', 'SS_M6', 'Deviance_M6',
'Alpha_M7', 'Delta_M7', 'Rho_M7', 'Beta_M7', 'SS_M7', 'Deviance_M7',
'Alpha_M8', 'Delta_M8', 'Rho_M8', 'Beta_M8', 'Epsilon_M8', 'SS_M8', 'Deviance_M8',
'Gamma_M9', 'SS_M9', 'Deviance_M9')
for (i in 2:ncol(altSubjectData)){altSubjectData[,i] = as.numeric(altSubjectData[,i])}
Compute Model Fit Index for Each Subject, for Each Alternative Model
altSubjectData$BIC_M1 = (altSubjectData$Deviance_M1) + log(65) * 4
altSubjectData$BIC_M2 = (altSubjectData$Deviance_M2) + log(65) * 4
altSubjectData$BIC_M3 = (altSubjectData$Deviance_M3) + log(65) * 4
altSubjectData$BIC_M4 = (altSubjectData$Deviance_M4) + log(65) * 5
altSubjectData$BIC_M5 = (altSubjectData$Deviance_M5) + log(65) * 5
altSubjectData$BIC_M6 = (altSubjectData$Deviance_M6) + log(65) * 5
altSubjectData$BIC_M7 = (altSubjectData$Deviance_M7) + log(65) * 4
altSubjectData$BIC_M8 = (altSubjectData$Deviance_M8) + log(65) * 5
altSubjectData$BIC_M9 = (altSubjectData$Deviance_M9) + log(65) * 1
Compare Model Performance
modelBIC = c(sum(subjectData$BIC), sum(altSubjectData$BIC_M1), sum(altSubjectData$BIC_M2), sum(altSubjectData$BIC_M3), sum(altSubjectData$BIC_M4), sum(altSubjectData$BIC_M5), sum(altSubjectData$BIC_M6), sum(altSubjectData$BIC_M7), sum(altSubjectData$BIC_M8), sum(altSubjectData$BIC_M9))
which(modelBIC == min(modelBIC))