Perceived-wealth-and-difficulty-modulation/Code/behavioralAnalyses.m

% Computes behavioral analyses on the dACC dataset.
% function developed by Demetrio Ferro. (C) CC BY-NC-SA 4.0

clear all; close all; clc;

savefolder='./Figures/';

dACC =load('../Data/dACC_all.mat');

Ntrs_dACC=arrayfun(@(cc)  dACC.data{cc}.numTrials, 1:length(dACC.data));

Nccs_dACC=length(dACC.data);
vars1_dACC=arrayfun(@(cc) dACC.data{cc}.vars(1,3), 1:length(dACC.data));
%unique_idx_dACC=find(diff([-1 Ntrs_dACC])~=0); % for consecutive sets with same number of trials, takes only the first
unique_idx_dACC=find(diff([-1 vars1_dACC])~=0); % for consecutive sets with same number vars1, takes only the first

dACC_s1=find(diff([-1 Ntrs_dACC(1:55)])~=0);
dACC_s2=find(diff([-1 Ntrs_dACC(56:129)])~=0);
dACC_s12=[dACC_s1 55+dACC_s2];

%% Loading main data from vars data structures
Vtl_dACC=struct2varsmat(dACC,5);
Vbl_dACC=struct2varsmat(dACC,6);
Vtr_dACC=struct2varsmat(dACC,7);
Vbr_dACC=struct2varsmat(dACC,8);
Ptl_dACC=struct2varsmat(dACC,12)./100;
Ptr_dACC=struct2varsmat(dACC,13)./100;
Ord_dACC=struct2varsmat(dACC,15);
Chlr_dACC=struct2varsmat(dACC,16);

%% Computing EV and Risk

EVl_dACC = Ptl_dACC.*Vtl_dACC + (1-Ptl_dACC).*Vbl_dACC;
EVr_dACC = Ptr_dACC.*Vtr_dACC + (1-Ptr_dACC).*Vbr_dACC;

EVf_dACC=nan(size(EVl_dACC));
EVs_dACC=nan(size(EVl_dACC));
EVf_dACC(Ord_dACC==1)=EVl_dACC(Ord_dACC==1);
EVf_dACC(Ord_dACC==2)=EVr_dACC(Ord_dACC==2);
EVs_dACC(Ord_dACC==1)=EVr_dACC(Ord_dACC==1);
EVs_dACC(Ord_dACC==2)=EVl_dACC(Ord_dACC==2);

Vtf_dACC=nan(size(Ord_dACC));
Vbf_dACC=nan(size(Ord_dACC));
Vts_dACC=nan(size(Ord_dACC));
Vbs_dACC=nan(size(Ord_dACC));
Vtf_dACC(Ord_dACC==1)=Vtl_dACC(Ord_dACC==1);
Vbf_dACC(Ord_dACC==1)=Vbl_dACC(Ord_dACC==1);
Vts_dACC(Ord_dACC==1)=Vtr_dACC(Ord_dACC==1);
Vbs_dACC(Ord_dACC==1)=Vbr_dACC(Ord_dACC==1);
Vts_dACC(Ord_dACC==2)=Vtl_dACC(Ord_dACC==2);
Vbs_dACC(Ord_dACC==2)=Vbl_dACC(Ord_dACC==2);
Vtf_dACC(Ord_dACC==2)=Vtr_dACC(Ord_dACC==2);
Vbf_dACC(Ord_dACC==2)=Vbr_dACC(Ord_dACC==2);

Chfs_dACC=nan(size(Ord_dACC));
Chfs_dACC(Ord_dACC==1)=Chlr_dACC(Ord_dACC==1);
Chfs_dACC(Ord_dACC==2)=(1-(Chlr_dACC(Ord_dACC==2)/2-1)*2);
ChC_dACC=nan(size(Ord_dACC));
ChC_dACC(EVf_dACC>=EVs_dACC)=double(Chfs_dACC(EVf_dACC>=EVs_dACC)==1);
ChC_dACC(EVf_dACC<=EVs_dACC)=double(Chfs_dACC(EVf_dACC<=EVs_dACC)==2);

Devlr_dACC=(EVl_dACC-EVr_dACC);
Devfs_dACC=(EVf_dACC-EVs_dACC);

Rl_dACC = Ptl_dACC.*(Vtl_dACC-EVl_dACC).^2 + (1-Ptl_dACC).*(Vbl_dACC-EVl_dACC).^2;
Rr_dACC = Ptr_dACC.*(Vtr_dACC-EVr_dACC).^2 + (1-Ptr_dACC).*(Vbr_dACC-EVr_dACC).^2;

Rf_dACC=nan(size(Ord_dACC));
Rs_dACC=nan(size(Ord_dACC));
Rf_dACC(Ord_dACC==1) = Rl_dACC(Ord_dACC==1);
Rf_dACC(Ord_dACC==2) = Rr_dACC(Ord_dACC==2);
Rs_dACC(Ord_dACC==1) = Rr_dACC(Ord_dACC==1);
Rs_dACC(Ord_dACC==2) = Rl_dACC(Ord_dACC==2);

%% Logit choice vs EV difference dACC sgACC

hff=figure('Units','Normalized','Position',[0 .33 .5 .33]);
subplot(1,2,1);
fitplotlogitglm((EVf_dACC(:)-EVs_dACC(:)),double((Chfs_dACC(:))==1),'r');
title('dACC: logit(p_{ch1}) = \beta_0+\beta_1 (EV_{1} - EV_{2})');
xlabel('EV_{1}-EV_{2}'); ylabel('p_{ch1} = probability of (ch=1)');

subplot(1,2,2);
fitplotlogitglm((EVr_dACC(:)-EVl_dACC(:)),double((Chlr_dACC(:))==2),'b');
title('dACC: logit(p_{chR}) = \beta_0+\beta_1 (EV_{R} - EV_{L})');
xlabel('EV_{R}-EV_{L}'); ylabel('p_{chR} = probability of (ch=R)');

saveas(hff,[savefolder 'Fig1_logit_ch_EV.png']);
saveas(hff,[savefolder 'Fig1_logit_ch_EV.svg']);


%% Logit choice vs Risk difference

hff=figure('Units','Normalized','Position',[0 0 .5 .33]);
subplot(1,2,1);
fitplotlogitglm((Rf_dACC(:)-Rs_dACC(:)),double((Chfs_dACC(:))==1),'r');
title('dACC: logit(p_{ch1})=\beta_0+\beta_1 (R_{1}-R_{2})');
xlabel('R_{1}-R_{2}'); ylabel('p_{ch1} = probability of (ch=1)');

subplot(1,2,2);
fitplotlogitglm((Rr_dACC(:)-Rl_dACC(:)),double((Chlr_dACC(:))==2),'b');
title('dACC: logit(p_{chR})=\beta_0+\beta_1 (R_{R}-R_{L})');
xlabel('R_{R}-R_{L}'); ylabel('p_{chR} = probability of (ch=R)');

saveas(hff,[savefolder '/Fig3_logit_ch_Risk.png']);
saveas(hff,[savefolder '/Fig3_logit_ch_Risk.svg']);


% ---------------------------------SUPPORT----FUNCTIONS----------------------------

%% Compute Logistic fits
function [gfit,ginfos]=fitplotlogitglm(X, Y, colorvs, alpha, linespec, linewidth, ...
    beta_offset, beta_color, plot_flag, Nbins, fontsize)
% plotmsem plots mean +/- sem with shaded, continuous error bars
% function developed by Demetrio Ferro. (C) CC BY-NC-SA 4.0
if nargin<11
    if nargin<10
        if nargin<9
            if nargin <7
                if nargin < 6
                    if nargin < 5
                        if nargin < 4
                            alpha=.1;
                        end
                        linespec='-';
                    end
                    linewidth=2;
                end
                beta_offset=[max(get(gca,'XLim')).*1.2, .1, .1];
            end
            plot_flag=1;
        end
        Nbins=50;
    end
    fontsize=6;
end
if nargin>=3 && ischar(colorvs) && contains('kbgcrmyw',colorvs(1))
    color=rem(floor((strfind('kbgcrmyw',colorvs(1))-1)*[.25 .5 1]),2);
    if length(colorvs)>1
        linespec=colorvs(2:end);
    end
elseif nargin>=3 && isa(colorvs,'double')
    color=colorvs(1:3);
else
    color=[0 0 0];
end
if nargin<8
    beta_color=color;
end

warning('off','all');

if length(Y)<=1
    szx=size(X);
    gfit.Coefficients.Estimate(1:(szx(szx>1)+1))=nan;
    ginfos.warn=1;
    ginfos.nobs=1;
    ginfos.warnmsg='One single observation is not enough for glm fit';

elseif length(Y)>1
    if size(X,2)>1 && size(X,1)>1 && plot_flag
        error('Plot function only works with one single regressor');
    elseif size(X,2)>1 && size(X,1)~=length(Y)
        X=X';
    end

    lastwarn('');
    gfit=fitglm(X,Y,'Distribution','binomial','link','logit');
    ginfos.nobs=sum(~isnan(Y));
    ginfos.warnmsg=lastwarn();
    if ~isempty(ginfos.warnmsg)
        ginfos.warn=1;
    else
        ginfos.warn=0;
    end

    if plot_flag
        gcci=gfit.coefCI;
        bp=gfit.Coefficients.pValue;
        bv=gfit.Coefficients.Estimate;
        bx = bv(1) + (X*bv(2));
        hx = 1 ./ (1 + exp(-sort(bx)));

        if Nbins>0
            nluqEV=linspace(floor(min(X)),ceil(max(X)),Nbins);
            [zq,xq]=nanmeanquantized(nluqEV,X,Y);

            plot(xq,zq,'.','color',color); hold on;
            plot(sort(X),hx,linespec,'color',color,'Linewidth',linewidth);
            plot([0 0],[0 1],'k-','LineWidth',1);
        end
        plot([nluqEV(1) 0],[1 1]./(1+exp(-bv(1))),'--','color',color);
        fill([xq flip(xq)],...
            [glmval(gcci(:,1),xq,'logit')' flip(glmval(gcci(:,2),xq,'logit')')],...
            color,'edgecolor','none','facealpha',alpha);
        %if length(nluqEV)>1; xlim(nluqEV([1 end])); daspect([nluqEV(end)-nluqEV(1) 1 1]); end
        text(beta_offset(1),beta_offset(2)+beta_offset(3)*2,...
            sprintf(['\\beta_0 = %1.3f^{' getpstars(bp(1)) '};'],bv(1)),'color',beta_color,'FontSize',fontsize);
        text(beta_offset(1),beta_offset(2)+beta_offset(3)*1,...
            sprintf(['\\beta_1 = %1.3f^{' getpstars(bp(2)) '};'],bv(2)),'color',beta_color,'FontSize',fontsize);
        text(beta_offset(1),beta_offset(2)+beta_offset(3)*0,...
            sprintf('n= %d;',length(X)),'color',beta_color,'FontSize',fontsize);
        box off
    end
end
end


function [ zq, xq ] = nanmeanquantized( x, y, z )
% find avg z at index y with value x
xq=nan(1,length(x)-1);
zq=nan(1,length(x)-1);

for ix=1:(length(x)-1)
    xq(ix)=x(ix)/2+x(ix+1)/2;
    if ix==(length(x)-1)
        zq(ix)=nanmean(z(y>=x(ix) & y<=x(ix+1)));
    else
        zq(ix)=nanmean(z(y>=x(ix) & y<x(ix+1)));
    end
end

end

function [ pstars ] = getpstars( pvalue )
if pvalue <0.001
    pstars=' ***';
elseif pvalue<0.01
    pstars=' ** ';
elseif pvalue<0.05
    pstars='  * ';
else
    pstars='  - ';
end
end

%% Get task variables
function ACCmat=struct2varsmat(ACCstruct, vars_index, start_offset,end_offset)
% The function takes as first input ACCstruct loaded from "../Data/dACC_all.mat" and provides task variables,
% respectively indicized by the second input vars_index, in trial range start_offset:end_offset
% function developed by Demetrio Ferro. (C) CC BY-NC-SA 4.0
if nargin<4
    if nargin<3
        start_offset=0;
    end
    end_offset=0;
end
Nccs=length(ACCstruct.data);
Ntrs=arrayfun(@(cc) ACCstruct.data{cc}.numTrials, 1:Nccs);

ACCmat=nan(Nccs,max(Ntrs));

for cc=1:Nccs
    Ntrcc=size(ACCstruct.data{cc}.vars,1);
    ACCmat(cc,:)=cat(1,ACCstruct.data{cc}.vars(:,vars_index),nan(max(Ntrs)-Ntrcc,1));
    ACCmat(cc,1:start_offset)=nan;
    ACCmat(cc,Ntrcc-end_offset+1:Ntrcc)=nan;
end
end