Estimating Free Parameters ************* Lesson ================ .. article-info:: :avatar: UCLA_Suit.jpg :avatar-link: https://www.decisionneurosciencelab.com/elijahgalvan :author: Elijah Galvan :date: September 1, 2023 :read-time: 11 min read :class-container: sd-p-2 sd-outline-muted sd-rounded-1 Goal During this Stage --------------- Now that we have our data, our first objective is to determine what the :bdg-success:`Free Parameters` are for each :bdg-success:`Subject`. How to Achieve this Goal ------------ .. dropdown:: Preallocating and Defining Functions .. tab-set:: .. tab-item:: Plain English 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. 1. :bdg-primary:`Trial` List - this should be in a certain order and have all of the trials that the subject has seen 2. :bdg-secondary:`Construct Value` Functions - this should relate :bdg-danger:`Decisions`, :bdg-primary:`Independent Variables`, and :bdg-primary:`Constants` to a number which encapsulates how much :bdg-danger:`Decisions` violate/follow a relevant norm 3. :bdg-warning:`Utility` Function - this should be a function of :bdg-secondary:`Construct Values` and :bdg-success:`Free Parameters` 4. Your Objective Function - a function which takes :bdg-success:`Free Parameters` and :bdg-danger:`Decisions` and returns the error between Expected :bdg-warning:`Utility` and Observed :bdg-warning:`Utility` 5. 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. 1. 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. 1. A trial-level data structure 2. A subject-level data structure .. tab-item:: R :: 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() .. tab-item:: MatLab :: 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(); .. tab-item:: Python :: 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() .. dropdown:: Define the :bdg-success:`Subject` Loop .. tab-set:: .. tab-item:: Plain English We're going to start our most superior ``for`` loop which iterates over :bdg-success:`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: 1. We need the data for this :bdg-success:`Subject`` 2. The :bdg-danger:`Decisions` for this :bdg-success:`Subject` need to be in the same order as the :bdg-primary:`Trial` List we use in our Objective Function 3. We need to determine what is going to be outputted .. dropdown:: So what are we starting with in this loop? A :bdg-success:`Subject` .. dropdown:: And what do we want to finish this loop with? :bdg-success:`Free Parameters` for this :bdg-success:`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. .. dropdown:: So what do we need to preallocate before this loop starts? An output for the :bdg-success:`Free Parameters` we'll estimate, along with any other subject information. Also, we'll output all trial-by-trial :bdg-success:`Subject` data that will be relevant later. Both of these are already done, nice. .. dropdown:: Then, what do we need to compute within this loop? :bdg-success:`Free Parameters` .. tab-item:: R :: 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 } .. tab-item:: MatLab :: 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 .. tab-item:: Python :: 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 .. dropdown:: Estimate :bdg-success:`Free Parameters` .. tab-set:: .. tab-item:: Plain English Now, we are going to answer the Determine Free Parameters demand placed on us in the :bdg-success:`Subject` loop, namely to estimate :bdg-success:`Free Parameters`. We first need to hand our Objective Function the :bdg-success:`Subject`'s data. Then, we need to store our data before we proceed to the next :bdg-success:`Subject`. .. dropdown:: So what are we starting with? :bdg-danger:`Decisions`, correctly ordered .. dropdown:: And what do we want to finish with? A single set of :bdg-success:`Free Parameters` .. dropdown:: So what do we need to preallocate? Nothing, we've already got everything we need. .. dropdown:: Then, what do we need to compute? We need to get data about what the model actually predicts based on the estimated :bdg-success:`Free Parameters`. .. tab-item:: R :: 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 } .. tab-item:: MatLab :: 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 .. tab-item:: Python :: 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 .. dropdown:: Determine Predicted :bdg-danger:`Decisions` for these :bdg-success:`Free Parameters` .. tab-set:: .. tab-item:: Plain English Now, we are going to answer the Determine Predictions demand placed on us. We have found the :bdg-success:`Subject`'s :bdg-success:`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! .. dropdown:: So what are we starting with? :bdg-success:`Free Parameters`, :bdg-danger:`Decisions`, and the :bdg-primary:`Trial` Set .. dropdown:: And what do we want to finish with? Predicted :bdg-danger:`Decisions` and the Model Error (which we will compute by comparing Predicted-and-Observed :bdg-danger:`Decisions`) A tip here, always name your columns immediately below your loop so that you don't forget what is what! .. dropdown:: So what do we need to preallocate? A vector for our predicted :bdg-danger:`Decisions`. .. dropdown:: Then, what do we need to compute? Nothing more. .. tab-item:: R :: 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') .. tab-item:: MatLab :: 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'}; .. tab-item:: Python :: 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 ---------------------- .. dropdown:: Preallocating and Defining Functions .. tab-set:: .. tab-item:: R :: 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 .. tab-item:: MatLab :: 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; .. tab-item:: Python :: 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 .. dropdown:: Define the :bdg-success:`Subject` Loop .. tab-set:: .. tab-item:: R :: for (i in 1:length(included_subjects)){ df = trialData[which(included_subjects[i] == trialData$Subject), ] #Estimate Free Parameters } .. tab-item:: MatLab :: for i = 1:length(included_subjects) df = trialData(trialData.Subject == included_subjects(i), :); % Estimate Free Parameters end .. tab-item:: Python :: 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 .. dropdown:: Estimate :bdg-success:`Free Parameters` .. tab-set:: .. tab-item:: R :: 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 } .. tab-item:: MatLab :: 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 .. tab-item:: Python :: 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 .. dropdown:: Determine Predicted :bdg-danger:`Decisions` for these :bdg-success:`Free Parameters` .. tab-set:: .. tab-item:: R :: 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') .. tab-item:: MatLab :: 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'}); .. tab-item:: Python :: 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 ------------------- .. dropdown:: Preallocating and Defining Functions .. tab-set:: .. tab-item:: R :: 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') .. tab-item:: MatLab :: 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'}; .. tab-item:: Python :: 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'] .. dropdown:: Define the :bdg-success:`Subject` and :bdg-primary:`Condition` Loops .. tab-set:: .. tab-item:: R :: 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 } } .. tab-item:: MatLab :: 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 .. tab-item:: Python :: 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 .. dropdown:: Estimate :bdg-success:`Free Parameters` for each :bdg-primary:`Condition` .. tab-set:: .. tab-item:: R :: 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 } } .. tab-item:: MatLab :: 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 .. tab-item:: Python :: 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 .. dropdown:: Determine Predicted :bdg-danger:`Decisions` for these :bdg-success:`Free Parameters` .. tab-set:: .. tab-item:: R :: 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') .. tab-item:: MatLab :: 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'}; .. tab-item:: Python :: 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 ------------------- .. dropdown:: Preallocating and Defining Functions .. tab-set:: .. tab-item:: R :: 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") }) } .. tab-item:: MatLab :: .. tab-item:: Python :: .. dropdown:: Estimate :bdg-success:`Free Parameters` and Determine Predicted :bdg-danger:`Decisions` .. tab-set:: .. tab-item:: R :: 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") .. tab-item:: MatLab :: .. tab-item:: Python :: Tutorial 4 - Li et al., 2022 ------------------- .. dropdown:: Preallocating and Defining Functions .. tab-set:: .. tab-item:: R :: 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) } .. tab-item:: MatLab :: .. tab-item:: Python :: .. dropdown:: Estimate :bdg-success:`Free Parameters` and Determine Predicted :bdg-danger:`Decisions` .. tab-set:: .. tab-item:: R :: 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") .. tab-item:: MatLab :: .. tab-item:: Python ::