function DataB=runRegressions(DataB,DesignMatrixID)

% This script runs the RT GLM's with the log RT's / distribution poisson link function log
%keyboard;
disp('Running GLMs')
for Isub=1:length(DataB.ID)
    
    
        %% What trials to include in GLM analysis
        switch DataB.GLMs(DesignMatrixID).dataType
            case 'SelfRTs'
                DataB.GLMs(DesignMatrixID).subject(Isub).IndexGLM           =ones(size(DataB.decision.self.reaction_time(:,Isub),1),1);
            case 'OtherRTs'
                DataB.GLMs(DesignMatrixID).subject(Isub).IndexGLM           =ones(size(DataB.decision.other.reaction_time(:,Isub),1),1);
            case 'SelfAccuracy'
                DataB.GLMs(DesignMatrixID).subject(Isub).IndexGLM          =ones(size(DataB.decision.self.correct(:,Isub),1),1);
            case 'OtherAccuracy'
                DataB.GLMs(DesignMatrixID).subject(Isub).IndexGLM          =ones(size(DataB.decision.other.correct(:,Isub),1),1);
            case 'BothChoices'
                DataB.GLMs(DesignMatrixID).subject(Isub).IndexGLM          =ones(size(DataB.decision.self.redchoice(:,Isub),1)+size(DataB.decision.other.redchoice(:,Isub),1),1);
            case 'SelfChoices'
                DataB.GLMs(DesignMatrixID).subject(Isub).IndexGLM          =ones(size(DataB.decision.self.redchoice(:,Isub),1),1);
            case 'OtherChoices'
                DataB.GLMs(DesignMatrixID).subject(Isub).IndexGLM          =ones(size(DataB.decision.other.redchoice(:,Isub),1),1);
            case 'Ratings'
                DataB.GLMs(DesignMatrixID).subject(Isub).IndexGLM           = ones(size(DataB.agent_rating.score(:,Isub),1),1); 
            case 'ChoiceDivergence'
                              DataB.GLMs(DesignMatrixID).subject(Isub).IndexGLM          =ones(size(DataB.decision.self.redchoice(:,Isub),1),1);  
                
        end
                IndexGLM=DataB.GLMs(DesignMatrixID).subject(Isub).IndexGLM==1;
                DataB.GLMs(DesignMatrixID).subject(Isub).IndexGLM=DataB.GLMs(DesignMatrixID).subject(Isub).IndexGLM==1;
        switch DesignMatrixID
            case 1                
                DataB.GLMs(DesignMatrixID).DesignMatrixLabels={'BlueEvidenceShared';'BlueEvidenceOccluded'};
                DesignMatrixraw=[DataB.evidence.Other(:,Isub) DataB.evidence.Overall(:,Isub)-DataB.evidence.Other(:,Isub)];
                DesignMatrix= DesignMatrixraw(IndexGLM,:); % Only have the to be included Trials in Design Matrix
            case 2
                DataB.GLMs(DesignMatrixID).DesignMatrixLabels={'OccludedAmountEvid';'OtherAmountEvid'};
                DesignMatrixraw=[(abs(DataB.evidence.Overall(:,Isub))-abs(DataB.evidence.Other(:,Isub))) abs(DataB.evidence.Other(:,Isub))];
                DesignMatrix= DesignMatrixraw(IndexGLM,:); % Only have the to be included Trials in Design Matrix
            case 3
                DataB.GLMs(DesignMatrixID).DesignMatrixLabels={'OccludedAmountEvid';'OtherAmountEvid'};
                DesignMatrixraw=[(abs(DataB.evidence.Overall(:,Isub))-abs(DataB.evidence.Other(:,Isub))) abs(DataB.evidence.Other(:,Isub))];
                DesignMatrix= DesignMatrixraw(IndexGLM,:); % Only have the to be included Trials in Design Matrix                
            case 4
                DataB.GLMs(DesignMatrixID).DesignMatrixLabels={'BlueEvidenceShared';'BlueEvidenceOccluded'};
                DesignMatrixraw=[DataB.evidence.Other(:,Isub) DataB.evidence.Overall(:,Isub)-DataB.evidence.Other(:,Isub)];
                DesignMatrix= DesignMatrixraw(IndexGLM,:); % Only have the to be included Trials in Design Matrix
            case 5
                DataB.GLMs(DesignMatrixID).DesignMatrixLabels={'OccludedAmountEvid';'OtherAmountEvid'};
                DesignMatrixraw=[(abs(DataB.evidence.Overall(:,Isub))-abs(DataB.evidence.Other(:,Isub))) abs(DataB.evidence.Other(:,Isub))];
                DesignMatrix= DesignMatrixraw(IndexGLM,:); % Only have the to be included Trials in Design Matrix
            case 6
                DataB.GLMs(DesignMatrixID).DesignMatrixLabels={'OccludedAmountEvid';'OtherAmountEvid'};
                DesignMatrixraw=[(abs(DataB.evidence.Overall(:,Isub))-abs(DataB.evidence.Other(:,Isub))) abs(DataB.evidence.Other(:,Isub))];
                DesignMatrix= DesignMatrixraw(IndexGLM,:); % Only have the to be included Trials in Design Matrix
            
            case 7   % Agent ability specific analyses
                %keyboard;
                AgentAbilityLvlis=nan(length(DataB.AgentAbilityLvl(:,Isub)),length(unique(DataB.AgentAbilityLvl(:))));
                ExistingAbilities=sort(unique(DataB.AgentAbilityLvl(:)),'ascend');
                for L=1:length(unique(DataB.AgentAbilityLvl(:)))
                AgentAbilityLvlis(DataB.AgentAbilityLvl(:,Isub)==ExistingAbilities(L),L)=1;
                end
                
                DataB.GLMs(DesignMatrixID).DesignMatrixLabels={'AgentAbilityLvl';'SharedBlueEvid_AL1';'OccludedBlueEvid_AL1';'SharedBlueEvid_AL2';'OccludedBlueEvid_AL2'};                
                DesignMatrixraw=[DataB.AgentAbilityLvl(:,Isub) [DataB.evidence.Other(:,Isub) DataB.evidence.Overall(:,Isub)-DataB.evidence.Other(:,Isub)].*repmat(AgentAbilityLvlis(:,1),1,2)...
                    [DataB.evidence.Other(:,Isub) DataB.evidence.Overall(:,Isub)-DataB.evidence.Other(:,Isub)].*repmat(AgentAbilityLvlis(:,2),1,2)];
                DesignMatrix= DesignMatrixraw(IndexGLM,:);
            case 8 % Sample specific analyses
            %keyboard;
            DataB.GLMs(DesignMatrixID).DesignMatrixLabels={'SharedSignal1';'SharedSignal2';'SharedSignal3';'SharedSignal4';'occludedSignal1';'occludedSignal2'...
                };   %;'STDSharedSignals';'STDoccludedSignals'
             DesignMatrixraw=[DataB.sharedSignals(:,:,Isub) DataB.occludedSignals(:,:,Isub)];
             DesignMatrixraw((DesignMatrixraw(:))==0)=nan;
             DesignMatrixraw=[DesignMatrixraw ];  %nanstd(DesignMatrixraw(:,1:4),0,2) nanstd(DesignMatrixraw(:,5:6),0,2)
                             DesignMatrix= DesignMatrixraw(IndexGLM,:);
%keyboard
            case 9 % absolute sample specific analyses
                DataB.GLMs(DesignMatrixID).DesignMatrixLabels={'SharedSignal1';'SharedSignal2';'SharedSignal3';'SharedSignal4';...
                    'occludedSignal1';'occludedSignal2'};
             %signFlip=2*(DataB.evidence.Overall(:,Isub)>0)-1;
             
             DesignMatrixraw=abs([DataB.sharedSignals(:,:,Isub) DataB.occludedSignals(:,:,Isub)]);             
             DesignMatrixraw((DesignMatrixraw(:)==0))=nan;
                               DesignMatrix= DesignMatrixraw(IndexGLM,:);
                              % keyboard;
                              
            case 10 
                
                AgentAbilityLvlis=nan(length(DataB.AgentAbilityLvl(:,Isub)),length(unique(DataB.AgentAbilityLvl(:))));
                ExistingAbilities=sort(unique(DataB.AgentAbilityLvl(:)),'ascend');
                for L=1:length(unique(DataB.AgentAbilityLvl(:)))
                AgentAbilityLvlis(DataB.AgentAbilityLvl(:,Isub)==ExistingAbilities(L),L)=1;
                end                
                
                DataB.GLMs(DesignMatrixID).DesignMatrixLabels={'AgentAbilityLvl';'SharedBlueEvid';'OccludedBlueEvid';'SharedBlueEvidByAbility';'OccludedBlueEvidByAbility'};                
                DesignMatrixraw=[DataB.AgentAbilityLvl(:,Isub) DataB.evidence.Other(:,Isub) DataB.evidence.Overall(:,Isub)-DataB.evidence.Other(:,Isub)...
                    [DataB.evidence.Other(:,Isub) DataB.evidence.Overall(:,Isub)-DataB.evidence.Other(:,Isub)].*normalise(DataB.AgentAbilityLvl(:,Isub))];
                DesignMatrix= DesignMatrixraw(IndexGLM,:);   
              case 11 % Sample specific analyses
            %keyboard;
            DataB.GLMs(DesignMatrixID).DesignMatrixLabels={'SharedSignal1';'SharedSignal2';'SharedSignal3';'SharedSignal4';'occludedSignal1';'occludedSignal2'...
                ;'STDSharedSignals';'STDoccludedSignals'};
             DesignMatrixraw=[DataB.sharedSignals(:,:,Isub) DataB.occludedSignals(:,:,Isub) ];
             DesignMatrixraw((DesignMatrixraw(:))==0)=nan;
             DesignMatrixraw=[DesignMatrixraw nanstd(DesignMatrixraw(:,1:4),0,2) nanstd(DesignMatrixraw(:,5:6),0,2)];
                             DesignMatrix= DesignMatrixraw(IndexGLM,:);

            case 12 % absolute sample specific analyses
                DataB.GLMs(DesignMatrixID).DesignMatrixLabels={'SharedSignal1';'SharedSignal2';'SharedSignal3';'SharedSignal4';'occludedSignal1';'occludedSignal2'...
                   };  % ;'STDSharedSignals';'STDoccludedSignals'
                DesignMatrixraw=[DataB.sharedSignals(:,:,Isub) DataB.occludedSignals(:,:,Isub) ];
                DesignMatrixraw((DesignMatrixraw(:))==0)=nan;
                DesignMatrixraw=[DesignMatrixraw ];  %nanstd(DesignMatrixraw(:,1:4),0,2) nanstd(DesignMatrixraw(:,5:6),0,2)
                DesignMatrix= DesignMatrixraw(IndexGLM,:);
                
                 case 13 
                     AgentAbilityLvlis=nan(length(DataB.AgentAbilityLvl(:,Isub)),length(unique(DataB.AgentAbilityLvl(:))));
                ExistingAbilities=sort(unique(DataB.AgentAbilityLvl(:)),'ascend');
                for L=1:length(unique(DataB.AgentAbilityLvl(:)))
                AgentAbilityLvlis(DataB.AgentAbilityLvl(:,Isub)==ExistingAbilities(L),L)=1;
                end         
                DataB.GLMs(DesignMatrixID).DesignMatrixLabels={'AgentLvl';'AL1_SharedSignal1';'AL1_SharedSignal2';'AL1_SharedSignal3';'AL1_SharedSignal4';...
                    'AL1_occludedSignal1';'AL1_occludedSignal2';'AL2_SharedSignal1';'AL2_SharedSignal2';'AL2_SharedSignal3';'AL2_SharedSignal4';...
                    'AL2_occludedSignal1';'AL2_occludedSignal2'};
                DesignMatrixraw=[DataB.AgentAbilityLvl(:,Isub)  [DataB.sharedSignals(:,:,Isub) DataB.occludedSignals(:,:,Isub)].*AgentAbilityLvlis(:,1) ...
                    [DataB.sharedSignals(:,:,Isub) DataB.occludedSignals(:,:,Isub)].*AgentAbilityLvlis(:,2)];
                DesignMatrixraw((DesignMatrixraw(:))==0)=nan;
                %DesignMatrixraw=[DesignMatrixraw nanstd(DesignMatrixraw(:,1:4),0,2) nanstd(DesignMatrixraw(:,5:6),0,2)];
                DesignMatrix= DesignMatrixraw(IndexGLM,:);
                
                case 14 
                    signFlip=2*(DataB.evidence.Other(:,Isub)>0)-1;
                DataB.GLMs(DesignMatrixID).DesignMatrixLabels={'SharedSignal1';'SharedSignal2';'SharedSignal3';'SharedSignal4';'counteroccludedSignal1';'counteroccludedSignal2'};
                DesignMatrixraw=[DataB.sharedSignals(:,:,Isub).*repmat(signFlip,1,4) (DataB.occludedSignals(:,:,Isub)).*repmat(-signFlip,1,2) ];
                DesignMatrix= DesignMatrixraw(IndexGLM,:);
               case 15
                                       signFlip=2*(DataB.evidence.Other(:,Isub)>0)-1;
                DataB.GLMs(DesignMatrixID).DesignMatrixLabels={'EvidenceShared';'counter EvidenceOccluded'};
                DesignMatrixraw=[abs(DataB.evidence.Other(:,Isub)) (DataB.evidence.Overall(:,Isub)-DataB.evidence.Other(:,Isub)).*-signFlip];
                DesignMatrix= DesignMatrixraw(IndexGLM,:);  
                
                 case 16
                     signFlip=2*(DataB.evidence.Other(:,Isub)>0)-1;
                     DataB.GLMs(DesignMatrixID).DesignMatrixLabels={'EvidenceShared';'counter EvidenceOccluded'};
                     DesignMatrixraw=[abs(DataB.evidence.Other(:,Isub)) (DataB.evidence.Overall(:,Isub)-DataB.evidence.Other(:,Isub)).*-signFlip];
                     DesignMatrix= DesignMatrixraw(IndexGLM,:);
                
                
            case 17
                DataB.GLMs(DesignMatrixID).DesignMatrixLabels={'BlueEvidenceShared';'BlueEvidenceOccluded';'SelfOther';'SOXBlueEvidS';'SOXBlueEvidO'};
                Evidences=repmat([DataB.evidence.Other(:,Isub) DataB.evidence.Overall(:,Isub)-DataB.evidence.Other(:,Isub)],2,1); % repeat to run same regressors for self and other choices
                Self_other_Dec_ID=[ones(length(DataB.evidence.Overall(:,Isub)),1);-1*ones(length(DataB.evidence.Overall(:,Isub)),1)];
                
                DesignMatrixraw=[Evidences Self_other_Dec_ID Evidences.*repmat(Self_other_Dec_ID,1,2)];
                DesignMatrix= DesignMatrixraw(IndexGLM,:); % Only have the to be included Trials in Design Matrix
                
                
            case 18
                DataB.GLMs(DesignMatrixID).DesignMatrixLabels={'BlueEvid_S Self';'BlueEvid_S Other';'BlueEvid_O_C Self';'BlueEvid_O_C Other'};
                
                DesignMatrixraw=[[DataB.evidence.Other(:,Isub);nan(length(DataB.evidence.Other(:,Isub)),1)] ...
                  [nan(length(DataB.evidence.Other(:,Isub)),1);DataB.evidence.Other(:,Isub) ]  ... 
                  [DataB.evidence.Overall(:,Isub)-DataB.evidence.Other(:,Isub);nan(length(DataB.evidence.Other(:,Isub)),1)] ...
                  [nan(length(DataB.evidence.Other(:,Isub)),1);DataB.evidence.Overall(:,Isub)-DataB.evidence.Other(:,Isub) ]];
                                    
                DesignMatrix= DesignMatrixraw(IndexGLM,:); % Only have the to be included Trials in Design Matrix
                
        end
                       
        % normalise design matrix by zscoring
        for NrC = 1:size(DesignMatrix,2)
            DesignMatrix(~isnan(DesignMatrix(:,NrC)),NrC)=zscore(DesignMatrix(~isnan(DesignMatrix(:,NrC)),NrC));
        end
        DesignMatrix(isnan(DesignMatrix))=0;   % After normalization, substitute NaN's with zeros
        DataB.GLMs(DesignMatrixID).normDesignMatrix(:,:,Isub)                                        =DesignMatrix;
        DataB.GLMs(DesignMatrixID).DesignMatrixCorrelations(:,:,Isub)                =corrcoef(DesignMatrix);
        
        % Determine Data type to be run in the GLM
        switch DataB.GLMs(DesignMatrixID).dataType
            case 'SelfRTs'
                y=DataB.decision.self.reaction_time(:,Isub) ;
            case 'OtherRTs'
                y=DataB.decision.other.reaction_time(:,Isub) ;
            case 'SelfChoices'
                y=1-DataB.decision.self.redchoice(:,Isub);
            case 'OtherChoices'
                y=1-DataB.decision.other.redchoice(:,Isub) ;
            case 'SelfAccuracy'
                y=DataB.decision.self.correct(:,Isub);
            case 'OtherAccuracy'
                y=DataB.decision.other.correct(:,Isub) ;
            case 'BothChoices'
                y=[1-DataB.decision.self.redchoice(:,Isub);1-DataB.decision.other.redchoice(:,Isub)];
                
            case 'BothRTs'
                %y=DataB.
            case 'Ratings'
                y=DataB.agent_rating.score(:,Isub) ;
            case 'ChoiceDivergence'
                y=DataB.decision.other.redchoice(:,Isub)~=DataB.decision.self.redchoice(:,Isub);
        end
        
        % Transform data manually, if necessary
        switch DataB.GLMs(DesignMatrixID).dataTransform
            case 'log'
                transy=log(y);
                y=transy;
                %disp(['log Transformation of ' DataB.GLMs(DesignMatrixID).dataType ' Data done'])
            case 'non'
                %disp(['No Transformation of ' DataB.GLMs(DesignMatrixID).dataType ' Data done'])
        end
        
        %% Run GLM , Allows for different link function, i.e. types of linear regressions.  This is so that RT data can be fit with gamma and choice data with binomial
        switch DataB.GLMs(DesignMatrixID).linkfunction
            case 'poisson'
                [betaweights(:,Isub),devf,stats{Isub}]=glmfit(DesignMatrix,y(IndexGLM),'poisson','link','log');
            case 'normal'
                [betaweights(:,Isub),devf,stats{Isub}]=glmfit(DesignMatrix,y(IndexGLM),'normal','link','identity');
            case 'binomial'
                [betaweights(:,Isub),devf,stats{Isub}]=glmfit(DesignMatrix,y(IndexGLM),'binomial','link','logit');
            case 'gamma'
                [betaweights(:,Isub),devf,stats{Isub}]=glmfit(DesignMatrix,y(IndexGLM),'gamma','link','log');
        end 
        
        %keyboard;
end
DataB.GLMs(DesignMatrixID).betaweights=betaweights;
DataB.GLMs(DesignMatrixID).stats=stats;
 


