Estimating Free Parameters
Lesson
Elijah Galvan
September 1, 2023
11 min read
Goal During this Stage
Now that we have our data, our first objective is to determine what the Free Parameters are for each Subject.
How to Achieve this Goal
Preallocating and Defining Functions
Before we start estimating free parameters, we need to check off a pretty simple list. There are several functions and variables that we should already have: in the event that you are using a seperate workspace from your simulation, make sure that these are included.
Trial List - this should be in a certain order and have all of the trials that the subject has seen
Construct Value Functions - this should relate Decisions, Independent Variables, and Constants to a number which encapsulates how much Decisions violate/follow a relevant norm
Utility Function - this should be a function of Construct Values and Free Parameters
Your Objective Function - a function which takes Free Parameters and Decisions and returns the error between Expected Utility and Observed Utility
Optimizer Inputs - the Initial Parameters and Boundaries for your Objective Function
Here, if we have a data set with subjects who are excluded, we need to know which subjects should be included in the analysis.
A list of all subject folders/files that tell us which data to retrieve for analysis - these must include the subject ID to be able to be identified
We also need to preallocate the output of our analysis as well as the raw trial-by-trial data.
A trial-level data structure
A subject-level data structure
all_subjects = read.csv2() #a file with all subject names
excluded_subjects = #
included_subjects = all_subjects$subjectID[-which(all_subjects$subjectID == excluded_subjects)]
parentfolder = '' #the parentfolder where the subject folder/file is
restoffilepath = '' #everything after the subject folder/file name
trialData = data.frame()
subjectData = data.frame()
all_subjects = readtable('all_subjects.csv'); % Assuming your file is named 'all_subjects.csv'
excluded_subjects = [1, 3, 5]; % Replace with the subject IDs to be excluded
included_subjects = all_subjects.subjectID(~ismember(all_subjects.subjectID, excluded_subjects));
parentfolder = ''; % Insert your parent folder path
restoffilepath = ''; % Insert the rest of the file path
trialData = table();
subjectData = table();
all_subjects = pd.read_csv('all_subjects.csv', sep=';') # Assuming your file is named 'all_subjects.csv'
excluded_subjects = [1, 3, 5] # Replace with the subject IDs to be excluded
included_subjects = all_subjects[~all_subjects['subjectID'].isin(excluded_subjects)]['subjectID']
parentfolder = '' # Insert your parent folder path
restoffilepath = '' # Insert the rest of the file path
trialData = pd.DataFrame()
subjectData = pd.DataFrame()
Define the Subject Loop
We’re going to start our most superior for loop which iterates over Subjects included in our analysis.
Before we start talking about estimating parameters, let’s just make sure that we have our ducks in a row:
We need the data for this Subject`
The Decisions for this Subject need to be in the same order as the Trial List we use in our Objective Function
We need to determine what is going to be outputted
So what are we starting with in this loop?
A Subject
And what do we want to finish this loop with?
Free Parameters for this Subject as well as all of the relevant variables for assessing our model. Namely, this would be either be the Negative Log-Likelihood or Sum of Squared Errors of our model predictions.
We also want to output any variables we think will be relevant for additional analyses at a trial-level or at a subject-level.
So what do we need to preallocate before this loop starts?
An output for the Free Parameters we’ll estimate, along with any other subject information. Also, we’ll output all trial-by-trial Subject data that will be relevant later.
Both of these are already done, nice.
Then, what do we need to compute within this loop?
Free Parameters
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)
#Determine Free Parameters
subjectData[i, ] = #to determine
trialData[start:end, ] = #to determine
}
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);
% Determine Free Parameters
subjectData(i, :) = %to determine
trialData(start:end, :) = %to determine
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)
# Determine Free Parameters
subjectData[i, :] = #to determine
trialData[start:end, :] = #to determine
Estimate Free Parameters
Now, we are going to answer the Determine Free Parameters demand placed on us in the Subject loop, namely to estimate Free Parameters. We first need to hand our Objective Function the Subject’s data. Then, we need to store our data before we proceed to the next Subject.
So what are we starting with?
Decisions, correctly ordered
And what do we want to finish with?
A single set of Free Parameters
So what do we need to preallocate?
Nothing, we’ve already got everything we need.
Then, what do we need to compute?
We need to get data about what the model actually predicts based on the estimated Free Parameters.
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)
#Just Added
result = fmincon(obj_function,x0 = initial_params, A = NULL, b = NULL, Aeq = NULL, beq = NULL,
lb = lower_bounds, ub = upper_bounds,
df = df)
# 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);
% Just Added
options = optimoptions('fmincon', 'Display', 'off');
result = fmincon(@(params) obj_function(params, df), initial_params, [], [], [], [], lower_bounds, upper_bounds, [], options);
% 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 = np.genfromtxt(datafile, delimiter=',', skip_header=1)
# Just Added
result = fmin_con(obj_function, x0=initial_params, bounds=(lower_bounds, upper_bounds), args = df)
# Determine Predictions
Determine Predicted Decisions for these Free Parameters
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!
So what are we starting with?
Free Parameters, Decisions, and the Trial Set
And what do we want to finish with?
Predicted Decisions and the Model Error (which we will compute by comparing Predicted-and-Observed Decisions)
A tip here, always name your columns immediately below your loop so that you don’t forget what is what!
So what do we need to preallocate?
A vector for our predicted Decisions.
Then, what do we need to compute?
Nothing more.
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)
result = fmincon(obj_function,x0 = initial_params, A = NULL, b = NULL, Aeq = NULL, beq = NULL,
lb = lower_bounds, ub = upper_bounds,
df = df)
#Just Added
closestPoint = which(as.numeric(freeParameters[,1]) == as.numeric(round(result$par[1])) & as.numeric(freeParameters[,2]) == as.numeric(round(result$par[2])))
df$Prediction = vector('numeric')
for (k in 1:length(df$Decisions)){
Utility = vector('numeric', length(Choices))
for (n in 1:length(Choices)){
Utility[n] = utility(parameter1 = results$par[1],
parameter2 = results$par[2],
construct1 = construct1(df$IV[k], df$Constant[k], Choices[n]),
construct2 = construct2(df$IV[k], df$Constant[k], Choices[n])),
construct3 = construct3(df$IV[k], df$Constant[k], Choices[n])
}
correct_choice = which(Utility == max(Utility))
if (length(correct_choice) > 1){
correct_choice = correct_choice[sample(correct_choice, 1)]
}
df$Prediction[k] = Choices[correct_choice]
}
model_NLL = -2 * log(sum(dnorm(df$Decision, mean = df$Prediction)))
model_SS = sum((df$Decision - df$Prediction)**2)
subjectData[i, ] = c(included_subjects[i], result$par[1], result$par[2], freeParameters$Strategy[closestPoint], model_NLL, model_SS)
#add any additional subject-level variables; if we have a priori clusters, you can include the strategy like we've done here
start = length(subjectData[, 1]) + 1
end = start + length(df$Decisions)
trialData[start:end, 1] = included_subjects[i]
trialData[start:end, 2] = df$IV
trialData[start:end, 3] = df$Constant
trialData[start:end, 4] = df$Decision
trialData[start:end, 5] = df$Prediction
}
colnames(subjectData) = c('SubjectID', 'Parameter1', 'Parameter2', 'Strategy', 'modelNLL', 'modelSS')
colnames(trialData) = c('SubjectID', 'IV', 'Constant', 'Decision', 'Prediction')
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);
result = fmincon(@obj_function, initial_params, [], [], [], [], lower_bounds, upper_bounds, df);
% Just Added
closestPoint = find(str2double(freeParameters(:,1)) == round(result(1)) & str2double(freeParameters(:,2)) == round(result(2)));
df.Prediction = zeros(size(df.Decisions));
for k = 1:length(df.Decisions)
Utility = zeros(size(Choices));
for n = 1:length(Choices)
Utility(n) = utility(result(1), result(2), construct1(df.IV(k), df.Constant(k), Choices(n)), construct2(df.IV(k), df.Constant(k), Choices(n)), construct3(df.IV(k), df.Constant(k), Choices(n)));
end
correct_choice = find(Utility == max(Utility));
if length(correct_choice) > 1
correct_choice = correct_choice(randi(length(correct_choice), 1));
end
df.Prediction(k) = Choices(correct_choice);
end
model_NLL = -2 * log(sum(normpdf(df.Decision, df.Prediction)));
model_SS = sum((df.Decision - df.Prediction).^2);
subjectData(i, :) = [included_subjects{i}, result(1), result(2), freeParameters.Strategy(closestPoint), model_NLL, model_SS];
start = size(subjectData, 1) + 1;
end_ = start + length(df.Decisions);
trialData(start:end_, 1) = included_subjects{i};
trialData(start:end_, 2) = df.IV;
trialData(start:end_, 3) = df.Constant;
trialData(start:end_, 4) = df.Decision;
trialData(start:end_, 5) = df.Prediction;
end
subjectData.Properties.VariableNames = {'SubjectID', 'Parameter1', 'Parameter2', 'Strategy', 'modelNLL', 'modelSS'};
trialData.Properties.VariableNames = {'SubjectID', 'IV', 'Constant', 'Decision', 'Prediction'};
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')
result = fmincon(obj_function, x0=initial_params, A=None, b=None, Aeq=None, beq=None, lb=lower_bounds, ub=upper_bounds, decisions=df['Decisions'])
# Just Added
closestPoint = np.where((freeParameters[:, 0].astype(float) == round(result[0])) & (freeParameters[:, 1].astype(float) == round(result[1])))[0]
df['Prediction'] = np.zeros(len(df['Decisions']))
for k in range(len(df['Decisions'])):
Utility = np.zeros(len(Choices))
for n in range(len(Choices)):
Utility[n] = utility(result[0], result[1], construct1(df['IV'][k], df['Constant'][k], Choices[n]), construct2(df['IV'][k], df['Constant'][k], Choices[n]), construct3(df['IV'][k], df['Constant'][k], Choices[n]))
correct_choice = np.where(Utility == max(Utility))[0]
if len(correct_choice) > 1:
correct_choice = np.random.choice(correct_choice, 1)
df['Prediction'][k] = Choices[correct_choice[0]]
model_NLL = -2 * np.log(np.sum(norm.pdf(df['Decision'], df['Prediction'])))
model_SS = np.sum((df['Decision'] - df['Prediction'])**2)
subjectData[i, :] = [included_subjects[i], result[0], result[1], freeParameters['Strategy'][closestPoint[0]], model_NLL, model_SS]
start = subjectData.shape[0] + 1
end_ = start + len(df['Decisions'])
trialData[start:end_, 0] = included_subjects[i]
trialData[start:end_, 1] = df['IV']
trialData[start:end_, 2] = df['Constant']
trialData[start:end_, 3] = df['Decision']
trialData[start:end_, 4] = df['Prediction']
subjectData.columns = ['SubjectID', 'Parameter1', 'Parameter2', 'Strategy', 'modelNLL', 'modelSS']
trialData.columns = ['SubjectID', 'IV', 'Constant', 'Decision', 'Prediction']
Tutorials
Tutorial 1 - van Baar, Chang, & Sanfey, 2019
Preallocating and Defining Functions
included_subjects = c(t(read.csv2('C:/Users/DELL/Downloads/tutorial1_Data/subjectsIncluded_batch1.csv', sep=',', header = F)),
t(read.csv2('C:/Users/DELL/Downloads/tutorial1_Data/subjectsIncluded_batch2.csv', sep=',', header = F)))
trialData = read.csv2("C:/Users/DELL/Downloads/tutorial1_Data/allDataLong.csv", sep =',')
trialData$Prediction = vector('numeric', length(trialData$Subject))
trialData$Strategy = vector('numeric', length(trialData$Subject))
trialData = trialData[-which(trialData$Investment == 0), ]
trialData = trialData[-which(is.na(as.numeric(trialData$Returned))),]
subjectData = data.frame()
trialList = read.csv2("C:/Users/DELL/Downloads/tutorial1_Data/trialSet.csv", sep =',', )
trialList = trialList[,2:3]
trialList$Believed_Multiplier = 4
trialList$Endowment = 10
subjects_batch1 = table2array(readtable('C:/Users/DELL/Downloads/tutorial1_Data/subjectsIncluded_batch1.csv', 'ReadVariableNames', false))';
subjects_batch2 = table2array(readtable('C:/Users/DELL/Downloads/tutorial1_Data/subjectsIncluded_batch2.csv', 'ReadVariableNames', false))';
included_subjects = [subjects_batch1, subjects_batch2];
trialData = readtable("C:/Users/DELL/Downloads/tutorial1_Data/allDataLong.csv", 'Delimiter', ',');
trialData.Prediction = zeros(height(trialData), 1);
trialData.Strategy = zeros(height(trialData), 1);
trialData = trialData(trialData.Investment ~= 0, :);
trialData = trialData(~isnan(str2double(trialData.Returned)), :);
subjectData = table();
trialList = readtable("C:/Users/DELL/Downloads/tutorial1_Data/trialSet.csv", 'Delimiter', ',');
trialList = trialList(:, 2:3);
trialList.Believed_Multiplier = 4;
trialList.Endowment = 10;
subjects_batch1 = pd.read_csv('C:/Users/DELL/Downloads/tutorial1_Data/subjectsIncluded_batch1.csv', header=None).transpose().values.flatten()
subjects_batch2 = pd.read_csv('C:/Users/DELL/Downloads/tutorial1_Data/subjectsIncluded_batch2.csv', header=None).transpose().values.flatten()
included_subjects = np.concatenate([subjects_batch1, subjects_batch2])
trialData = pd.read_csv("C:/Users/DELL/Downloads/tutorial1_Data/allDataLong.csv", sep=',')
trialData['Prediction'] = np.zeros(len(trialData['Subject']))
trialData['Strategy'] = np.zeros(len(trialData['Subject']))
trialData = trialData[trialData['Investment'] != 0]
trialData = trialData.dropna(subset=['Returned'])
subjectData = pd.DataFrame()
trialList = pd.read_csv("C:/Users/DELL/Downloads/tutorial1_Data/trialSet.csv", sep=',')
trialList = trialList.iloc[:, 1:3]
trialList['Believed_Multiplier'] = 4
trialList['Endowment'] = 10
Define the Subject Loop
for (i in 1:length(included_subjects)){
df = trialData[which(included_subjects[i] == trialData$Subject), ]
#Estimate Free Parameters
}
for i = 1:length(included_subjects)
df = trialData(trialData.Subject == included_subjects(i), :);
% Estimate Free Parameters
end
for i in range(len(included_subjects)):
df = trialData[trialData['Subject'] == included_subjects[i]]
result = minimize(obj_function, x0=initial_params, bounds=list(zip(lower_bounds, upper_bounds)), args=(df,))
#Estimate Free Parameters
Estimate Free Parameters
for (i in 1:length(included_subjects)){
df = trialData[which(included_subjects[i] == trialData$Subject), ]
result = fmincon(obj_function,x0 = initial_params, A = NULL, b = NULL, Aeq = NULL, beq = NULL,
lb = lower_bounds, ub = upper_bounds,
df = df)
# Determine Predictions
}
for i = 1:length(included_subjects)
df = trialData(trialData.Subject == included_subjects(i), :);
result = fmincon(@obj_function, initial_params, [], [], [], [], lower_bounds, upper_bounds, [], df);
% Determine Predictions
end
for i in range(len(included_subjects)):
df = trialData[trialData['Subject'] == included_subjects[i]]
result = minimize(obj_function, x0=initial_params, bounds=list(zip(lower_bounds, upper_bounds)), args=(df,))
# Determine Predictions
Determine Predicted Decisions for these Free Parameters
for (i in 1:length(included_subjects)){
df = trialData[which(included_subjects[i] == trialData$Subject), ]
result = fmincon(obj_function,x0 = initial_params, A = NULL, b = NULL, Aeq = NULL, beq = NULL,
lb = lower_bounds, ub = upper_bounds,
df = df)
closestPoint = which(as.numeric(freeParameters[,1], 3) == (round(result$par[1]/2, 2))*2 &
((round((as.numeric(freeParameters[,2]) + 0.1)*5, 2))/5) - 0.1 == ((round((result$par[2] + 0.1)*5, 2))/5) - 0.1)
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)}
Utility = vector('numeric', length(Choices))
for (n in 1:length(Choices)){
Utility[n] = utility(theta = result$par[1],
phi = result$par[2],
guilt = guilt(I, B, Choices[n], M),
inequity = inequity(I, M, Choices[n], E),
payout = payout_maximization(I, M, Choices[n]))
}
correct_choice = which(Utility == max(Utility))
if (length(correct_choice) > 1){
correct_choice = correct_choice[sample(1:length(correct_choice), 1)]
}
df$Prediction[k] = Choices[correct_choice]
}
model_NLL = -2 * log(sum(dnorm(as.numeric(df$Returned), mean = df$Prediction)))
model_SS = sum((as.numeric(df$Returned) - df$Prediction)**2)
subjectData[i, 1:6] = c(included_subjects[i], result$par[1], result$par[2], freeParameters$Strategy[closestPoint], model_NLL, model_SS)
trialData$Prediction[which(included_subjects[i] == trialData$Subject)] = df$Prediction
trialData$Strategy[which(included_subjects[i] == trialData$Subject)] = freeParameters$Strategy[closestPoint]
}
colnames(subjectData) = c('SubjectID', 'Theta', 'Phi', 'Strategy', 'modelNLL', 'modelSS')
subjectData = zeros(length(included_subjects), 6);
for i = 1:length(included_subjects)
subj_id = included_subjects(i);
df = trialData(trialData.Subject == subj_id, :);
result = fmincon(@obj_function, initial_params, [], [], [], [], lower_bounds, upper_bounds, [], df);
closestPoint = find(round(freeParameters(:,1), 3) == round(result(1)/2, 2)*2 & ...
(round((freeParameters(:,2) + 0.1)*5, 2)/5 - 0.1 == round((result(2) + 0.1)*5, 2)/5 - 0.1), 1);
for k = 1:length(df.Returned)
I = df.Investment(k);
M = df.Multiplier(k);
R = df.Returned(k);
B = 4;
E = 10;
if I > 10
Choices = 0:(I*M)/10:(I*M);
else
Choices = 0:(I * M);
end
Utility = zeros(size(Choices));
for n = 1:length(Choices)
Utility(n) = utility(result(1), result(2), ...
guilt(I, B, Choices(n), M), ...
inequity(I, M, Choices(n), E), ...
payout_maximization(I, M, Choices(n)));
end
[~, correct_choice] = max(Utility);
df.Prediction(k) = Choices(correct_choice);
end
model_NLL = -2 * log(sum(normpdf(df.Returned, df.Prediction)));
model_SS = sum((df.Returned - df.Prediction).^2);
subjectData(i, :) = [subj_id, result(1), result(2), freeParameters.Strategy(closestPoint), model_NLL, model_SS];
trialData.Prediction(trialData.Subject == subj_id) = df.Prediction;
trialData.Strategy(trialData.Subject == subj_id) = freeParameters.Strategy(closestPoint);
end
subjectData = array2table(subjectData, 'VariableNames', {'SubjectID', 'Theta', 'Phi', 'Strategy', 'modelNLL', 'modelSS'});
subjectData = np.zeros((len(included_subjects), 6))
for i, subj_id in enumerate(included_subjects):
df = trialData[trialData['Subject'] == subj_id]
result = minimize(obj_function, x0=initial_params, bounds=list(zip(lower_bounds, upper_bounds)), args=(df,))
closestPoint = np.where((np.round(freeParameters.iloc[:, 0], 3) == np.round(result.x[0]/2, 2)*2) &
(np.round((freeParameters.iloc[:, 1] + 0.1)*5, 2)/5 - 0.1 == np.round((result.x[1] + 0.1)*5, 2)/5 - 0.1))[0]
for k in range(len(df['Returned'])):
I = df['Investment'].iloc[k]
M = df['Multiplier'].iloc[k]
R = df['Returned'].iloc[k]
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(theta=result.x[0],
phi=result.x[1],
guilt=guilt(I, B, Choices[n], M),
inequity=inequity(I, M, Choices[n], E),
payout=payout_maximization(I, M, Choices[n]))
correct_choice = np.argmax(Utility)
df['Prediction'].iloc[k] = Choices[correct_choice]
model_NLL = -2 * np.log(np.sum(norm.pdf(df['Returned'], loc=df['Prediction'])))
model_SS = np.sum((df['Returned'] - df['Prediction'])**2)
subjectData[i, :] = [subj_id, result.x[0], result.x[1], freeParameters['Strategy'].iloc[closestPoint].values[0], model_NLL, model_SS]
trialData.loc[trialData['Subject'] == subj_id, 'Prediction'] = df['Prediction']
trialData.loc[trialData['Subject'] == subj_id, 'Strategy'] = freeParameters['Strategy'].iloc[closestPoint].values[0]
subjectData = pd.DataFrame(subjectData, columns=['SubjectID', 'Theta', 'Phi', 'Strategy', 'modelNLL', 'modelSS'])
Tutorial 2 - Galvan & Sanfey, 2024
Preallocating and Defining Functions
all_subjects = read.csv2('C:/Users/DELL/Downloads/prolific_export_648c19e5420c9b10a79589a4.csv', sep = ',')
potentially_excluded_subjects = read.csv2('C:/Users/DELL/Downloads/prolific_export_6437e471bec005411b1503ea.csv', sep = ',')
k = 0
for (i in 1:length(potentially_excluded_subjects$Participant.id)){
k = k + sum(potentially_excluded_subjects$Participant.id[i] == all_subjects$Participant.id)
}
actual_excluded_subjects = vector('character', k)
j = 0
for (i in 1:length(potentially_excluded_subjects$Participant.id)){
if (sum(potentially_excluded_subjects$Participant.id[i] == all_subjects$Participant.id) == 1) {
j = j + 1
actual_excluded_subjects[j] = which(potentially_excluded_subjects$Participant.id[i] == all_subjects)
}
}
included_subjects = all_subjects[-actual_excluded_subjects]
parentfolder = 'C:/Users/DELL/rdb_all_modded/' #the parentfolder where the subject folder/file is
restoffilepath = '.csv' #everything after the subject folder/file name
trialData = data.frame()
subjectData = data.frame()
conditions = c('merit', 'entitlement', 'corruption', 'luck')
all_subjects = readtable('C:/Users/DELL/Downloads/prolific_export_648c19e5420c9b10a79589a4.csv', 'Delimiter', ',');
potentially_excluded_subjects = readtable('C:/Users/DELL/Downloads/prolific_export_6437e471bec005411b1503ea.csv', 'Delimiter', ',');
k = 0;
for i = 1:length(potentially_excluded_subjects.Participant_id)
k = k + sum(potentially_excluded_subjects.Participant_id(i) == all_subjects.Participant_id);
end
actual_excluded_subjects = cell(1, k);
j = 0;
for i = 1:length(potentially_excluded_subjects.Participant_id)
if sum(potentially_excluded_subjects.Participant_id(i) == all_subjects.Participant_id) == 1
j = j + 1;
actual_excluded_subjects{j} = find(potentially_excluded_subjects.Participant_id(i) == all_subjects.Participant_id);
end
end
included_subjects = all_subjects(~ismember(1:length(all_subjects.Participant_id), cell2mat(actual_excluded_subjects)), :);
parentfolder = 'C:/Users/DELL/rdb_all_modded/';
restoffilepath = '.csv';
trialData = table();
subjectData = table();
conditions = {'merit', 'entitlement', 'corruption', 'luck'};
all_subjects = pd.read_csv('C:/Users/DELL/Downloads/prolific_export_648c19e5420c9b10a79589a4.csv', delimiter=',')
potentially_excluded_subjects = pd.read_csv('C:/Users/DELL/Downloads/prolific_export_6437e471bec005411b1503ea.csv', delimiter=',')
k = 0
for i in range(len(potentially_excluded_subjects['Participant.id'])):
k += sum(potentially_excluded_subjects['Participant.id'][i] == all_subjects['Participant.id'])
actual_excluded_subjects = [None] * k
j = 0
for i in range(len(potentially_excluded_subjects['Participant.id'])):
if sum(potentially_excluded_subjects['Participant.id'][i] == all_subjects['Participant.id']) == 1:
j += 1
actual_excluded_subjects[j-1] = potentially_excluded_subjects['Participant.id'][i]
included_subjects = all_subjects[~all_subjects['Participant.id'].isin(actual_excluded_subjects)]
parentfolder = 'C:/Users/DELL/rdb_all_modded/'
restoffilepath = '.csv'
trialData = pd.DataFrame()
subjectData = pd.DataFrame()
conditions = ['merit', 'entitlement', 'corruption', 'luck']
Define the Subject and Condition Loops
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)
thetaPerCondition = vector('numeric', length(conditions))
phiPerCondition = vector('numeric', length(conditions))
strategyPerCondition = vector('numeric', length(conditions))
SSPerCondition = vector('numeric', length(conditions))
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
#Determine Free Parameters for Each Condition
}
}
for i = 1:length(included_subjects)
datafile = fullfile(parentfolder, [included_subjects{i}, restoffilepath]);
fullData = readtable(datafile);
thetaPerCondition = zeros(1, length(conditions));
phiPerCondition = zeros(1, length(conditions));
strategyPerCondition = zeros(1, length(conditions));
SSPerCondition = zeros(1, length(conditions));
for j = 1:length(conditions)
df = fullData(fullData.condition == conditions(j), [49, 40:48, 33]);
df.redistributionRate = df.redistributionRate / 100;
% Determine Free Parameters for Each Condition
end
end
for i in range(len(included_subjects)):
datafile = os.path.join(parentfolder, f"{included_subjects[i]}{restoffilepath}")
fullData = pd.read_csv(datafile)
thetaPerCondition = np.zeros(len(conditions))
phiPerCondition = np.zeros(len(conditions))
strategyPerCondition = np.zeros(len(conditions))
SSPerCondition = np.zeros(len(conditions))
for j in range(len(conditions)):
df = fullData[fullData['condition'] == conditions[j]][['49', '40', '41', '42', '43', '44', '45', '46', '47', '48', '33']]
df['redistributionRate'] = df['redistributionRate'] / 100
# Determine Free Parameters for Each Condition
end
Estimate Free Parameters for each 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)
thetaPerCondition = vector('numeric', length(conditions))
phiPerCondition = vector('numeric', length(conditions))
strategyPerCondition = vector('numeric', length(conditions))
SSPerCondition = vector('numeric', length(conditions))
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
#Just Added
result = fmincon(obj_function,x0 = initial_params, A = NULL, b = NULL, Aeq = NULL, beq = NULL,
lb = lower_bounds, ub = upper_bounds,
df = df)
thetaPerCondition[j] = result$par[1]
phiPerCondition[j] = result$par[2]
closestPoint = which(as.numeric(freeParameters[,1]) == round(result$par[1], 2) & as.numeric(freeParameters[,2]) == round(result$par[2], 2))
strategyPerCondition[j] = freeParameters$Strategy[closestPoint]
#Determine Predictions
}
}
for i = 1:length(included_subjects)
datafile = fullfile(parentfolder, [included_subjects{i}, restoffilepath]);
fullData = readtable(datafile);
thetaPerCondition = zeros(1, length(conditions));
phiPerCondition = zeros(1, length(conditions));
strategyPerCondition = zeros(1, length(conditions));
SSPerCondition = zeros(1, length(conditions));
for j = 1:length(conditions)
df = fullData(fullData.condition == conditions(j), [49, 40:48, 33]);
df.redistributionRate = df.redistributionRate / 100;
% Just Added
result = fmincon(@(params)obj_function(params, df), initial_params, [], [], [], [], lower_bounds, upper_bounds);
thetaPerCondition(j) = result(1);
phiPerCondition(j) = result(2);
closestPoint = find(freeParameters.theta == round(result(1), 2) & freeParameters.phi == round(result(2), 2));
strategyPerCondition(j) = freeParameters.Strategy(closestPoint);
% Determine Predictions
end
end
for i in range(len(included_subjects)):
datafile = os.path.join(parentfolder, f"{included_subjects[i]}{restoffilepath}")
fullData = pd.read_csv(datafile)
thetaPerCondition = np.zeros(len(conditions))
phiPerCondition = np.zeros(len(conditions))
strategyPerCondition = np.zeros(len(conditions))
SSPerCondition = np.zeros(len(conditions))
for j in range(len(conditions)):
df = fullData[fullData['condition'] == conditions[j]][[48] + list(range(39, 4)) + [32]]
df['redistributionRate'] = df['redistributionRate'] / 100
# Just Added
result = fmincon(obj_function, initial_params, df);
thetaPerCondition[j] = result[0]
phiPerCondition[j] = result[1]
closestPoint = np.where((freeParameters['theta'] == round(result[0], 2)) & (freeParameters['phi'] == round(result[1], 2)))[0]
strategyPerCondition[j] = freeParameters['Strategy'][closestPoint].values[0]
# Determine Predictions
Determine Predicted Decisions for these Free Parameters
for (i in 1:length(included_subjects)){
datafile = paste(parentfolder, included_subjects[i], restoffilepath, sep = '') # produces a character vector 'parentfolder/included_subjects[i]restoffilepath'
fullData = read.csv2(datafile)
thetaPerCondition = vector('numeric', length(conditions))
phiPerCondition = vector('numeric', length(conditions))
strategyPerCondition = vector('numeric', length(conditions))
SSPerCondition = vector('numeric', length(conditions))
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
result = fmincon(obj_function,x0 = initial_params, A = NULL, b = NULL, Aeq = NULL, beq = NULL,
lb = lower_bounds, ub = upper_bounds,
df = df)
thetaPerCondition[j] = result$par[1]
phiPerCondition[j] = result$par[2]
closestPoint = which(as.numeric(freeParameters[,1]) == round(result$par[1], 2) & as.numeric(freeParameters[,2]) == round(result$par[2], 2))
strategyPerCondition[j] = freeParameters$Strategy[closestPoint]
#Just Added
df$predictedRR = vector('numeric')
df$predictedOutcome = vector('numeric')
df$actualOutcome = vector('numeric')
for (k in 1:length(df$redistributionRate)){
Utility = vector('numeric', length(Choices))
for (n in 1:length(Choices)){
Utility[n] = utility(theta = thetaPerCondition[j],
phi = phiPerCondition[j],
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]),
Payout = payout(new_value(df[k, 1], choices[n]), df[k, 1], choices[n]))
}
correct_choice = which(Utility == max(Utility))
df$predictedRR[k] = Choices[correct_choice[sample(length(correct_choice), 1)]]
df$predictedOutcome[k] = new_value(df$myself[k], df$predictedRR[k])
df$actualOutcome[k] = new_value(df$myself[k], df$redistributionRate[k])
}
SSPerCondition[j] = sum((df$predictedOutcome - df$actualOutcome)**2)
fullData[which(fullData$condition == conditions[j]), 76] = df$predictedOutcome
fullData[which(fullData$condition == conditions[j]), 77] = df$predictedRR
fullData[which(fullData$condition == conditions[j]), 78] = df$actualOutcome
fullData[which(fullData$condition == conditions[j]), 79] = freeParameters$Strategy[closestPoint]
}
subjectData[i, 1:17] = c(included_subjects[i], thetaPerCondition, phiPerCondition, strategyPerCondition, SSPerCondition)
start = length(subjectData[, 1]) + 1
end = start + length(fullData$redistributionRate)
trialData[start:end, 1] = included_subjects[i]
trialData[start:end, 2] = fullData$redistributionRate/100
trialData[start:end, 3] = fullData[,78]
trialData[start:end, 4] = fullData[,77]
trialData[start:end, 5] = fullData[,76]
trialData[start:end, 6] = fullData$myself
trialData[start:end, 7] = fullData[, 79]
trialData[start:end, 9] = rep(seq(1, 4), each = length(df$redistributionRate))
trialData[start:end, 10] = seq(1, length(fullData$redistributionRate))
trialData[start:end, 11] = fullData$condition
trialData[start:end, 12:20] = fullData[, c40:48]
}
colnames(subjectData) = c('SubjectID', 'thetaMerit', 'thetaEntitlement', 'thetaCorruption', 'thetaLuck', 'phiMerit', 'phiEntitlement', 'phiCorruption', 'phiLuck', 'strategyMerit', 'strategyEntitlement', 'strategyCorruption', 'strategyLuck', 'SSMerit', 'SSEntitlement', 'SSCorruption', 'SSLuck')
colnames(trialData) = c('SubjectID', 'observedTaxRate', 'observedOutcome', 'predictedTaxRate', 'predictedOutcome', 'initialAllocation', 'strategy', 'blockNumber', 'trialNumber', 'condition',
'P1', 'P2', 'P3', 'P4', 'P5', 'P6', 'P7', 'P8', 'P9')
for i = 1:length(included_subjects)
datafile = fullfile(parentfolder, [included_subjects{i}, restoffilepath]);
fullData = readtable(datafile);
thetaPerCondition = zeros(1, length(conditions));
phiPerCondition = zeros(1, length(conditions));
strategyPerCondition = zeros(1, length(conditions));
SSPerCondition = zeros(1, length(conditions));
for j = 1:length(conditions)
df = fullData(fullData.condition == conditions(j), [49, 40:48, 33]);
df.redistributionRate = df.redistributionRate / 100;
result = fmincon(@(params)obj_function(params, df), initial_params, [], [], [], [], lower_bounds, upper_bounds);
thetaPerCondition(j) = result(1);
phiPerCondition(j) = result(2);
closestPoint = find(freeParameters.theta == round(result(1), 2) & freeParameters.phi == round(result(2), 2));
strategyPerCondition(j) = freeParameters.Strategy(closestPoint);
df.predictedRR = zeros(length(df.redistributionRate), 1);
df.predictedOutcome = zeros(length(df.redistributionRate), 1);
df.actualOutcome = zeros(length(df.redistributionRate), 1);
for k = 1:length(df.redistributionRate)
Utility = zeros(length(Choices), 1);
for n = 1:length(Choices)
Utility(n) = utility(thetaPerCondition(j), phiPerCondition(j), ...
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)), ...
payout(new_value(df(k, 1), Choices(n)), df(k, 1), Choices(n)));
end
[~, correct_choice] = max(Utility);
df.predictedRR(k) = Choices(correct_choice);
df.predictedOutcome(k) = new_value(df.myself(k), df.predictedRR(k));
df.actualOutcome(k) = new_value(df.myself(k), df.redistributionRate(k));
end
SSPerCondition(j) = sum((df.predictedOutcome - df.actualOutcome).^2);
fullData(fullData.condition == conditions(j), 76) = df.predictedOutcome;
fullData(fullData.condition == conditions(j), 77) = df.predictedRR;
fullData(fullData.condition == conditions(j), 78) = df.actualOutcome;
fullData(fullData.condition == conditions(j), 79) = freeParameters.Strategy(closestPoint);
end
subjectData(i, 1:17) = [included_subjects{i}, thetaPerCondition, phiPerCondition, strategyPerCondition, SSPerCondition];
start = size(subjectData, 1) + 1;
stop = start + length(fullData.redistributionRate) - 1;
trialData(start:stop, 1) = included_subjects{i};
trialData(start:stop, 2) = fullData.redistributionRate / 100;
trialData(start:stop, 3) = fullData(:, 78);
trialData(start:stop, 4) = fullData(:, 77);
trialData(start:stop, 5) = fullData(:, 76);
trialData(start:stop, 6) = fullData.myself;
trialData(start:stop, 7) = fullData(:, 79);
trialData(start:stop, 9) = repelem(1:4, length(df.redistributionRate));
trialData(start:stop, 10) = 1:length(fullData.redistributionRate);
trialData(start:stop, 11) = fullData.condition;
trialData(start:stop, 12:20) = fullData(:, 40:48);
end
subjectData.Properties.VariableNames = {'SubjectID', 'thetaMerit', 'thetaEntitlement', 'thetaCorruption', 'thetaLuck', ...
'phiMerit', 'phiEntitlement', 'phiCorruption', 'phiLuck', 'strategyMerit', 'strategyEntitlement', ...
'strategyCorruption', 'strategyLuck', 'SSMerit', 'SSEntitlement', 'SSCorruption', 'SSLuck'};
trialData.Properties.VariableNames = {'SubjectID', 'observedTaxRate', 'observedOutcome', 'predictedTaxRate', ...
'predictedOutcome', 'initialAllocation', 'strategy', 'blockNumber', 'trialNumber', 'condition', ...
'P1', 'P2', 'P3', 'P4', 'P5', 'P6', 'P7', 'P8', 'P9'};
for i in range(len(included_subjects)):
datafile = os.path.join(parentfolder, f"{included_subjects[i]}{restoffilepath}")
fullData = pd.read_csv(datafile)
thetaPerCondition = np.zeros(len(conditions))
phiPerCondition = np.zeros(len(conditions))
strategyPerCondition = np.zeros(len(conditions))
SSPerCondition = np.zeros(len(conditions))
for j in range(len(conditions)):
df = fullData[fullData['condition'] == conditions[j]][[48] + list(range(39, 4)) + [32]]
df['redistributionRate'] = df['redistributionRate'] / 100
result = fmincon(obj_function, initial_params, df);
thetaPerCondition[j] = result[0]
phiPerCondition[j] = result[1]
closestPoint = np.where((freeParameters['theta'] == round(result[0], 2)) & (freeParameters['phi'] == round(result[1], 2)))[0]
strategyPerCondition[j] = freeParameters['Strategy'][closestPoint].values[0]
df['predictedRR'] = np.zeros(len(df['redistributionRate']))
df['predictedOutcome'] = np.zeros(len(df['redistributionRate']))
df['actualOutcome'] = np.zeros(len(df['redistributionRate']))
for k in range(len(df['redistributionRate'])):
Utility = np.zeros(len(Choices))
for n in range(len(Choices)):
Utility[n] = utility(thetaPerCondition[j], phiPerCondition[j],
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]),
payout(new_value(df.iloc[k, 0], Choices[n]), df.iloc[k, 0], Choices[n]))
correct_choice = np.argmax(Utility)
df['predictedRR'].iloc[k] = Choices[correct_choice]
df['predictedOutcome'].iloc[k] = new_value(df['myself'].iloc[k], df['predictedRR'].iloc[k])
df['actualOutcome'].iloc[k] = new_value(df['myself'].iloc[k], df['redistributionRate'].iloc[k])
SSPerCondition[j] = np.sum((df['predictedOutcome'] - df['actualOutcome'])**2)
fullData.loc[fullData['condition'] == conditions[j], 75] = df['predictedOutcome']
fullData.loc[fullData['condition'] == conditions[j], 76] = df['predictedRR']
fullData.loc[fullData['condition'] == conditions[j], 77] = df['actualOutcome']
fullData.loc[fullData['condition'] == conditions[j], 78] = freeParameters['Strategy'][closestPoint].values[0]
subjectData.iloc[i, 0:17] = [included_subjects[i], *thetaPerCondition, *phiPerCondition, *strategyPerCondition, *SSPerCondition]
start = len(subjectData) + 1
stop = start + len(fullData['redistributionRate'])
trialData.loc[start:stop, 'SubjectID'] = included_subjects[i]
trialData.loc[start:stop, 'observedTaxRate'] = fullData['redistributionRate'] / 100
trialData.loc[start:stop, 'observedOutcome'] = fullData.iloc[:, 77]
trialData.loc[start:stop, 'predictedTaxRate'] = fullData.iloc[:, 75]
trialData.loc[start:stop, 'predictedOutcome'] = fullData.iloc[:, 76]
trialData.loc[start:stop, 'initialAllocation'] = fullData['myself']
trialData.loc[start:stop, 'strategy'] = fullData.iloc[:, 78]
trialData.loc[start:stop, 'blockNumber'] = np.tile(np.arange(1, 5), len(df['redistributionRate']))
trialData.loc[start:stop, 'trialNumber'] = np.arange(1, len(fullData['redistributionRate']) + 1)
trialData.loc[start:stop, 'condition'] = fullData['condition']
trialData.loc[start:stop, 'P1':'P9'] = fullData.iloc[:, 39:47].values
Tutorial 3 - Crockett et al., 2014
Preallocating and Defining Functions
library(R.matlab) #data are stored in .mat files so we'll need this to work while in R
parentfolder = 'C:/Users/DELL/Downloads/fMRI_choice_data/run' #the parentfolder where the subject folder/file is
included_subjects = list.files(paste(parentfolder, '1', sep = ''), full.names = F) #get subject files in run 1
included_subjects = gsub('.mat', '', gsub('BC_1', '', included_subjects)) #remove run 1 identifier and .mat suffix
subjectData = data.frame()
grab_data = function(subject){
if (sum(grepl(pattern = subject, x = list.files(paste(parentfolder, '2/', sep = ''), full.names = F))) == 0){
return("skipping")
}
for (r in 1:2){
datafile = paste(parentfolder, r, '/BC_', r, subject, '.mat', sep = '')
temp = readMat(datafile)
if (r == 1){
choices = temp["choice"]
df = data.frame(temp["trialinfo"])
df = cbind(df[1:76, ], as.data.frame(choices)[1:76, ])
} else {
choices = temp["choice"]
temp = data.frame(temp["trialinfo"])
temp = cbind(temp[1:76, ], as.data.frame(choices)[1:76, ])
df = rbind(df, temp) #some have duplicate values
}
}
df[which(is.na(df[,11])), 11] = sample(c(1,2), size = sum(is.na(df[,11])), replace = T)
df[, 10] = (df[, 10] * -1) + 2 # 1 is harm 0 is help now
df[, 11] = (df[, 11] * -1) + 2 # 1 is self 0 is other now
colnames(df) = c("trialNumber", "moneyLess", "moneyDifference", "moneyMore", "shocksLess", "shocksDifference", "shocksMore", "isLeft", "runNumber", "isSelf", "choseHarm")
return(df)
}
generate_predictions = function(df, result){
df$ProbHarm = 0
Epsilon = result$par[5]
Gamma = result$par[6]
for (k in 1:length(df[,1])){
shocksMore = df$shocksMore[k]
shocksLess = df$shocksLess[k]
moneyMore = df$moneyMore[k]
moneyLess = df$moneyLess[k]
isLeft = as.logical(df$isLeft[k])
isSelf = as.logical(df$isSelf[k])
if (isSelf) {Kappa = result$par[1]; Beta = result$par[3]} else {Kappa = result$par[2]; Beta = result$par[4]}
Utility = utility(Harm = harm(shocksMore, shocksLess),
Payout = payout(moneyMore, moneyLess),
kappa = Kappa)
df$ProbHarm[k] = max(min(probability(Beta, Epsilon, Gamma, Utility, isLeft), 0.9999999999), 0.00000000001)
}
return(df)
}
obj_function = 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 = params[6]
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)))
}
}
#also need to include the correct number of parameter constraints
initial_params = c(0.5, 0.5, 4, 4, 0.25, 0)
lower_bounds = c(0, 0, 0, 0, 0, -0.5)
upper_bounds = c(1, 1, 100, 100, 0.5, 0.5)
#and overcome some errors that seem to be limited to using Log-Odds with fmincon
optimize = function(df){
tryCatch({
fmincon(obj_function, x0 = initial_params, lb = lower_bounds, ub = upper_bounds, df = df, optimMethod = "MLE")
}, error = function(e){
fmincon(obj_function, x0 = initial_params, lb = lower_bounds, ub = upper_bounds, df = df, optimMethod = "OLS")
})
}
Estimate Free Parameters and Determine Predicted Decisions
for (i in 1:length(included_subjects)){
df = grab_data(included_subjects[i])
if (is.character(df)){next}
result = optimize(df)
# Just Added
df = generate_predictions(df, result)
model_SS = sum((df$choseHarm - df$ProbHarm)**2)
model_NLL = -2*sum(df$choseHarm * log(df$ProbHarm) + (1 - df$choseHarm) * log(1 - df$ProbHarm))
subjectData[i, 1:9] = c(included_subjects[i], result$par[1], result$par[2], result$par[3], result$par[4], result$par[5], result$par[6],
model_SS, model_NLL)
df = cbind(data.frame(SubjectID = rep(included_subjects[i], length(df$trialNumber))), df)
if (i == 1) {trialData = df} else {trialData = rbind(trialData, df)}
}
colnames(subjectData) = c("subjectID", "Kappa_Self", "Kappa_Other", "Beta_Self", "Beta_Other", "Epsilon", "Gamma", "SS", "Deviance")
Tutorial 4 - Li et al., 2022
Preallocating and Defining Functions
trialData = read.csv2("C:/Users/DELL/Downloads/Data/Data/HPP_fMRI_beh_data_for_lmm.csv", sep = ',')
trialData = trialData[which(trialData$trail_type==3), c(1,6,7,13,14,21,22,36)]
colnames(trialData) = c('SubjectID', 'a0', 'b0', 'a1', 'b1', 'a2', 'b2', 'Chose1')
trialData$Chose1 = trialData$Chose1 - 1
trialData$Prob1 = 0
included_subjects = levels(factor(trialData$SubjectID))
subjectData = data.frame()
grab_data = function(subject){
eturn(trialData[which(trialData$SubjectID == subject), -1])
}
addPredictions = function(trialData, subject, predictions){
trialData$Prob1[which(trialData$SubjectID == subject)] = predictions
return(trialData)
}
Estimate Free Parameters and Determine Predicted Decisions
for (i in 1:length(included_subjects)){
df = grab_data(included_subjects[i])
result = optimize(obj_function, initial_params, lower_bounds, upper_bounds, df)
# Just Added
df$Prob1 = generatePredictions(result$par, 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))
subjectData[i, 1:9] = c(included_subjects[i], result$par[1], result$par[2], result$par[3], result$par[4], result$par[5], result$par[6],
model_SS, model_NLL)
trialData = addPredictions(trialData, included_subjects[i], df$Prob1)
}
colnames(subjectData) = c("subjectID", "Alpha", "Delta", "Rho", "Beta", "Epsilon", "Gamma", "SS", "Deviance")