% Study 1 clc, clear clear global rng(0,'twister'); global soft_out options = optimset;% optimization specs options = optimset(options,'Display', 'off','PlotFcns', {@optimplotx @optimplotfval}, ... 'Algorithm', 'active-set', 'FinDiffType', 'central'); options.PlotFcns = []; D=dir('*.csv');% array of subject files to loop through subjStay=[]; subjProbe=[]; AIC=[]; for s=1:length({D.name}) %% Data organizing format longg opts = detectImportOptions(D(s).name); opts = setvartype(opts, 'condition', 'char'); % Load subject-specific data dat = readtable(D(s).name, opts); % Color-to-condition assignment % Shape-to-condition assignment dat.Properties.VariableNames(15)={'shapes'}; shapesCB = erase(split(dat.shapes(1))', ["[", "]", "'", ","]); choice=dat.choice(~isnan(dat.choice)); trialType=dat.ProbeType(~isnan(dat.choice)); rewLev=dat.rewardLevel(~isnan(dat.choice)); cons=dat.payoffCons(~isnan(dat.choice));% consensus level on each trial chosUnchosOPT=zeros(height(choice),3);% array to plug in chosen & unchosen trial option + outcome condition conds={'1', '2','3', '4', '1_partial', '2_partial'};% consensus & win/lose conditions for cb=1:6 chosUnchosOPT(strcmp(dat.imageLeft, shapesCB{cb}),1)=cb;% add IDs of available shape options chosUnchosOPT(strcmp(dat.imageRight, shapesCB{cb}),2)=cb; chosUnchosOPT(strcmp(dat.condition, conds(cb)),3)=cb+6; end colprob={'42', '24', '13', '31'};% probe color options for cb=1:4 chosUnchosOPT(strcmp(dat.condition, colprob(cb)),1)=str2double(colprob{cb}(1))+6; chosUnchosOPT(strcmp(dat.condition, colprob(cb)),2)=str2double(colprob{cb}(2))+6; end chosUnchosREW=[dat.rewardLeft(~isnan(dat.choice)), dat.rewardRight(~isnan(dat.choice))];% payoff on trial for each option % Convert to Chosen & Unchosen options and rewards for i = 1:height(choice) if choice(i)==39 chosUnchosOPT(i,1:2)=[chosUnchosOPT(i,2),chosUnchosOPT(i,1)]; chosUnchosREW(i,1:2)=[chosUnchosREW(i,2),chosUnchosREW(i,1)]; end if chosUnchosREW(i,1)-1)'; postRating = dat.Var19(dat.Var19>-1)'; behavStay=bigStay./numTrials; for c=1:8 choice_probe(c)=length(find(chosUnchosOPT(:,1)==c & trialType(:,1)~=0))/length(find((chosUnchosOPT(:,1)==c | chosUnchosOPT(:,2)==c) & trialType(:,1)~=0)); rewProbe(c)=nanmean(chosUnchosREW(chosUnchosOPT(:, 1)==c,1)); end subjSpecs{s}={trialType, chosUnchosOPT, chosUnchosREW, preRating, cons, consWinLev};% Store subject-specific variable arrays % Cell array organizing data for plotting means, participants are rows, measures are columns subjStay(s,:)=nanmean(behavStay);% Stay probabilities in each condition, order is CW, CL, DW, DL, PW, PL, for high and low reward distributions subjProbe(s,:)=[nanmean(choice_probe(:, 1:2)), nanmean(choice_probe(:, 4:5)), choice_probe(:,3), ... choice_probe(:,6:8)];% Probe choice props for shapes by rew and freq levels, (HH, LH, HL, LL), followed by consensus pref totRew(s,:)=[rewStay(1,:), rewStay(2,:), [mean(rewProbe(1:2)),mean(rewProbe(:,4:5)), mean(rewProbe(:,3)), mean(rewProbe(:,6))]];% Rewards received in each Stay & probe condition totTrials(s,:)=[numTrials(1,:), numTrials(2,:)];% Trials in each Stay condition subjRat{s,5}=postRating-preRating; %% Fit & Compare Models k=[4,3,3,2];% # of parameters for AIC for mods=1:4 paramRec=0; numPar=repmat([0.5; 0.01; 0.99], 1, k(mods)-2); [params, nll, ~, ~]=fminsearchbnd(@(fp) fitparams(fp, subjSpecs{s}, mods, paramRec),[numPar(1,:),0.5,5],[numPar(2,:),0.01,1],[numPar(3,:), 0.99,10],options); subjParams{s,mods}=params; AIC(s,mods)=2*k(mods)+2*(nll); end end save behavSpecs subj* AIC %% Parameter Recovery paramRec=input('Enter 1 for Parameter Recovery: '); iter=1000; if paramRec == 1 % Randomly generated parameter values qlr=(0.99-0.01).*rand(iter,1) + 0.01;% learning rate for money ulr=(0.99-0.01).*rand(iter,1) + 0.01;% learning rate for action imitation vlr=(0.99-0.01).*rand(iter,1) + 0.01;% learning rate for social reward dn=(10-1).*rand(iter,1) + 1; % Decision noice theta=[qlr, ulr, vlr, dn]; modParam=[]; m=[1 5 9 13]; for pr=1:iter % generate experiment events using randomly drawn participant r=randi([1,length({D.name})]); specs=subjSpecs{r}; rep=[1 2 3]; for x=1:length(rep)% increase training trials specs{rep(x)}=repmat(specs{rep(x)}, 10,1); end specs{5}=randi(2, 1500,1)+randi(2, 1500,1)-2;% generate choices by 2 other agents for mods=1:length(m) % simulate weights & decisions given parameters & model fitparams(theta(pr,:),specs, mods, paramRec); simChoice{pr,mods}=soft_out; % recover parameters paramRec=-1; [params, nll, ~, ~]=fminsearchbnd(@(fp) fitparams(fp, simChoice{pr,mods}, mods, paramRec),[0.3,0.3,0.3,5],[0.01,0.01,0.01,1],[0.99,0.99,0.99,10],options); modParam(pr, m(mods):m(mods)+3)=params; paramRec=1; end end for mods=1:length(m) h=modParam(:,m(mods):m(mods)+3); for pr=1:width(theta) [r,p]=corr(h(:, pr), theta(:,pr)); paramCorrR(mods,pr)=r; paramCorrP(mods,pr)=p; end end paramRecovCorr=array2table(paramCorrP, 'RowNames', {'hybrid', 'imitation', 'social', 'baseline'},'VariableNames', {'$lr', 'imitationlr', 'sociallr', 'decision noise'}); save paramRecov simChoice paramCorr* theta modParam paramRecovCorr %% Model Predictions probeTrials=[]; stayProbs=[]; mPt=[1 9 17 25]; mSp=[1 7 13 19]; for x=1:length(mPt) for y=1:height(simChoice) % Probe trials hold=simChoice{y,x}; hold=hold{2}(hold{1}~=0, 1:2); pT=zeros(2,8); for z=1:8 pT(1,z)=size(find(hold(:,1)==z),1); pT(2,z)=size(find(hold(:,1:2)==z),1); end hold=simChoice{y,x}; hold=hold{2}(hold{1}==0, 1:3); hold(:,3)=hold(:,3)-6; stayTrials=zeros(2,6); for i=1:height(hold)-1 if hold(i+1,1) == hold(i,1) || hold(i+1,2) == hold(i,1) % if chosen option is repeated on the next trial stayTrials(1, hold(i,3))=stayTrials(1, hold(i,3))+1;% counter if hold(i+1,1) == hold(i,1) %if stayed stayTrials(2, hold(i,3))=stayTrials(2, hold(i,3))+1;% counter end end end probeTrials(y,mPt(x):mPt(x)+7)=pT(1,:)./pT(2,:); stayProbs(y,mSp(x):mSp(x)+5)=stayTrials(2,:)./stayTrials(1,:); end end end %% Models function nll=fitparams(fp, specs, mods, paramRec) global soft_out % specs: 1. trialType, 2. chosUnchosOPT, 3. chosUnchosREW, 4. rating, 5. cons, 6. consWinLev freqVals = ones(1,12); if mods<3 % hybrid and imitation actVals = zeros(1,12); payVals = zeros(1,12)+0.5; sumVals=actVals+payVals; else sumVals = zeros(1,12)+0.5; end % social and monetary initVals=sumVals; M=zeros(height(specs{1})); modelprob=zeros(size(M,1),1); specsO=specs; for i = 1:height(M) M(i,1)=sumVals(specs{2}(i,1))+sqrt(2*log(i)/freqVals(specs{2}(i,1))); M(i,2)=sumVals(specs{2}(i,2))+sqrt(2*log(i)/freqVals(specs{2}(i,2))); modelprob(i)=exp(fp(length(fp))*M(i,1))/(exp(fp(length(fp))*M(i,1))+exp(fp(length(fp))*M(i,2)));% softmax with noise parameter if paramRec==1 chosResp = randsample([1 0], 1, true, [modelprob(i) 1-modelprob(i)]);% generate response based on model probs if chosResp==0% chosen and unchosen response contingencies specs{2}(i,1)=specsO{2}(i,2); specs{3}(i,1)=specsO{3}(i,2); specs{2}(i,2)=specsO{2}(i,1); specs{3}(i,2)=specsO{3}(i,1); specs{5}(i)=2-specsO{5}(i); end if specs{5}(i)==2; specs{2}(i,3)=7; elseif specs{5}(i)==0; specs{2}(i,3)=9; elseif specs{5}(i)==1; specs{2}(i,3)=11; end if specs{3}(i,1)