phase analysis script

pjkohler · Jun 26, 2018 · 610085e · 610085e
1 parent 523a21a
commit 610085e
Showing 1 changed file with 337 additions and 0 deletions.
diff --git a/phase_analysis_m2D3D.m b/phase_analysis_m2D3D.m
@@ -0,0 +1,337 @@
+%% SCRIPT FOR DOING PHASE ANALYSIS
+
+codeFolder = '/Users/kohler/code';
+rcaCodePath = sprintf('%s/git/rcaBase',codeFolder);
+addpath(genpath(rcaCodePath));
+addpath(genpath(sprintf('%s/git/mrC',codeFolder)));
+addpath(genpath(sprintf('%s/git/schlegel/matlab_lib',codeFolder)));
+addpath(genpath('~/code/git/export_fig/'))
+setenv('DYLD_LIBRARY_PATH','')
+
+close all; clear superRCA;
+
+%% PREPARE
+
+% we now hardcode the paths to the two data sets
+dateStr = datestr(clock,26);
+dateStr(strfind(dateStr,'/')) ='';
+
+% make subject list
+topFolder = '/Volumes/Denali_4D2/kohler/EEG_EXP/DATA/motion2D3D/';
+
+expList = {'exp1','exp2','exp3','exp4'};
+condsToUse = [1,3,5,7]; % use full reference conditions
+
+idList = [];
+for e = 1:length(expList);
+    expFolder = subfolders(sprintf('%s/*%s*',topFolder,expList{e}),1);
+    if e == 1
+        subjList = subfolders(sprintf('%s/*nl-*',expFolder{1}),1);
+    else
+    end
+    tempList = subfolders(sprintf('%s/*nl-*',expFolder{1}),1);
+    for s = 1:length(tempList);
+        [~,curSub]=fileparts(tempList{s});
+        idIdx = strfind(curSub,'nl-');
+        curSub = curSub(idIdx:(idIdx+6));
+        if e == 1
+            idList = [idList;{curSub}];
+        else
+            if any(~cell2mat(cellfun(@(x) isempty(strfind(x,curSub)),subjList,'uni',false)))
+                continue;
+            else
+                subjList = [subjList;tempList(s)];
+                idList = [idList;{curSub}];
+            end
+        end
+    end 
+end
+
+% go deeper into the folder names
+for f = 1:length(subjList)
+    tempFolders = subfolders(subjList{f},1);
+    tempFolders = tempFolders(cellfun(@(x) isempty(strfind(x,'mff')),tempFolders));
+    folderNames{f} = sprintf('%s/Exp_TEXT_HCN_128_Avg',tempFolders{end});
+end
+
+saveFilePath = '/Volumes/Denali_4D2/kohler/EEG_EXP/DATA/motion2D3D/figures/phase_analysis';
+
+binsToUse=1:10; % indices of bins to include in analysis (the values must be present in the bin column of all DFT/RLS exports)
+freqsToUse= [2,4]; % indices of frequencies to include in analysis (the values must be present in the frequency column of all DFT/RLS exports)
+trialsToUse = []; %1:10; % subset of trials to use for analysis (if set to false or empty, all trials will be used)
+nReg=7; % RCA regularization constant (7-9 are typical values, but see within-trial eigenvalue plot in rca output)
+nComp=5; % number of RCs that you want to look at (3-5 are good values, but see across-trial eigenvalue plot in rca output)
+chanToCompare = 75; % channel to use for a performance evaluation, can be []
+dataType = 'RLS'; % can also be 'DFT' if you have DFT exports
+rcPlotStyle = 'matchMaxSignsToRc1'; % not req'd. see 'help rcaRun', can be: 'matchMaxSignsToRc1' (default) or 'orig'
+forceSourceData = false;
+trialError = false;
+
+%% RC analysis 
+condNames = {'HorizontalDisparity', 'VerticalDisparity'};
+for f = 1:length(freqsToUse)
+    superRCA(f) = rcaSweep(folderNames,binsToUse,freqsToUse(f),condsToUse,trialsToUse,nReg,nComp,dataType,chanToCompare,1,rcPlotStyle,forceSourceData);
+    export_fig(sprintf('%s/SuperSetHoriVert_%0.0f_cov.pdf',saveFilePath,freqsToUse(f)),'-pdf','-transparent','-opengl',gcf);
+    close gcf;
+end
+zeroData = cellfun(@(x) any(x(:)==0), cat(1,superRCA(:).data));
+
+if any(zeroData(:))
+    error('data values exactly zero, this should not happen');
+else
+end
+
+%% COMPUTE VALUES FOR PLOTTING
+
+keepConditions = true;
+errorType = 'SEM';
+doNR = false(2,6,4); % 2 freqs, 5 RCs, 2 conditions
+doNR(:,1,:) = true; % do fitting for first RC, second and fourth harmonic, all conditions
+
+for f = 1:length(superRCA)
+        fprintf('\n ... rc no. %0.0f ...\n',f);
+        rcStruct = aggregateData(superRCA(f),keepConditions,errorType,trialError,doNR);
+        % RC
+        superRCA(f).stats.Amp = squeeze(rcStruct.ampBins);
+        superRCA(f).stats.SubjectAmp = squeeze(rcStruct.subjectAmp);
+        superRCA(f).stats.ErrLB = squeeze(rcStruct.ampErrBins(:,:,:,:,1));
+        superRCA(f).stats.ErrUB = squeeze(rcStruct.ampErrBins(:,:,:,:,2));
+        superRCA(f).stats.NoiseAmp = squeeze(rcStruct.ampNoiseBins);
+        superRCA(f).stats.SubjectNoiseAmp = squeeze(rcStruct.subjectAmpNoise);
+        % Naka-Rushton
+        superRCA(f).stats.NR_Params = squeeze(rcStruct.NakaRushton.Params);
+        superRCA(f).stats.NR_R2 = squeeze(rcStruct.NakaRushton.R2);
+        superRCA(f).stats.NR_JKSE = squeeze(rcStruct.NakaRushton.JackKnife.SE);
+        superRCA(f).stats.NR_JKParams = squeeze(rcStruct.NakaRushton.JackKnife.Params);
+        superRCA(f).stats.hModel = rcStruct.NakaRushton.hModel;
+        % t-values
+        superRCA(f).stats.tSqrdP = squeeze(rcStruct.tSqrdP);
+        superRCA(f).stats.tSqrdSig = squeeze(rcStruct.tSqrdSig);
+        superRCA(f).stats.tSqrdVal = squeeze(rcStruct.tSqrdVal);
+end
+% shut down parallel pool, which was used for fitting Naka-Rushton
+delete(gcp('nocreate'));
+clc;
+
+ellipsePlot(superRCA(1))
+
+%% MAKE FIGURE
+close all
+rcNum = 1;
+lWidth = 1.5;
+fSize = 12;
+gcaOpts = {'tickdir','out','ticklength',[0.0500,0],'box','off','fontsize',fSize,'fontname','Helvetica','linewidth',lWidth};
+cBrewer = load('colorBrewer.mat');
+mainColors = [cBrewer.rgb20(5,:); cBrewer.rgb20(7,:)];
+
+topoVals = [superRCA(1).A(:,rcNum);superRCA(2).A(:,rcNum)];
+rcaColorBar = [min(topoVals),max(topoVals)];
+newExtreme = round(max(abs(rcaColorBar(:,f)))*5)./5;
+%rcaColorBar = [-newExtreme,newExtreme*1.001];
+rcaColorBar = [-0.4,0.4001];
+plotLabel = {'A','B','C','D'};
+
+for f=1:length(freqsToUse)
+    binVals = cellfun(@(x) str2num(x), superRCA(f).settings.binLabels);
+    logStep = diff(reallog(binVals(1:2))); % step size
+    xMin = reallog(binVals(1))-logStep*.5;
+    xMax = reallog(binVals(end))+logStep*2.5; % add 2.5 steps
+    extraBins = arrayfun(@(x) exp(reallog(binVals(end))+x), [logStep,logStep*2]);
+
+    % make Naka-Rushton values
+    NRmodel = superRCA(f).stats.hModel;
+    NRset = squeeze(superRCA(f).stats.NR_Params(:,rcNum,:));          
+    NRerrs = superRCA(f).stats.NR_JKSE;
+    NRmodel = superRCA(f).stats.hModel;
+
+    % compute new signal values, averaged over bins
+    [rcaDataReal,rcaDataImag] = getRealImag(superRCA(f).data);
+    rcaDataReal = cellfun(@(x) squeeze(nanmean(x(:,rcNum,:),3)),rcaDataReal,'uni',false);
+    rcaDataReal = cell2mat(permute(rcaDataReal,[3,2,1]));
+    rcaDataImag = cellfun(@(x) squeeze(nanmean(x(:,rcNum,:),3)),rcaDataImag,'uni',false);
+    rcaDataImag = cell2mat(permute(rcaDataImag,[3,2,1]));
+    realBinMean = squeeze(nanmean(rcaDataReal));
+    imagBinMean = squeeze(nanmean(rcaDataImag));
+    % compute new noise values, averaged over bins
+    [noiseLoReal,noiseLoImag] = getRealImag(superRCA(f).noiseData.lowerSideBand);
+    [noiseHiReal,noiseHiImag] = getRealImag(superRCA(f).noiseData.lowerSideBand);
+    noiseReal = cellfun(@(x,y) (x+y)./2, noiseLoReal,noiseHiReal, 'uni',false);
+    noiseImag = cellfun(@(x,y) (x+y)./2, noiseLoImag,noiseHiImag, 'uni',false);
+    noiseReal = cellfun(@(x) squeeze(nanmean(x(:,rcNum,:),3)),noiseReal,'uni',false);
+    noiseReal = cell2mat(permute(noiseReal,[3,2,1]));
+    noiseImag = cellfun(@(x) squeeze(nanmean(x(:,rcNum,:),3)),noiseImag,'uni',false);
+    noiseImag = cell2mat(permute(noiseImag,[3,2,1]));
+    realBinMeanNoise = squeeze(nanmean(noiseReal));
+    imagBinMeanNoise = squeeze(nanmean(noiseImag));
+
+    % compute projected vector mean
+    % move subjects to first dim
+    realVector = permute(rcaDataReal,[2,1,3]);
+    imagVector = permute(rcaDataImag,[2,1,3]);
+    project_amps = vectorProjection(realVector,imagVector);
+    % compute mean of projected vector amplitude, over bins
+    project_amps(:,end+1,:) = nanmean(project_amps,2);
+    % if all values are NaN
+    if any(sum(squeeze(all(isnan(project_amps),2)),2)>0)
+        nan_subs = find(sum(squeeze(all(isnan(project_amps),2)),2)>0);
+        for z = 1:length(nan_subs)
+            msg = ...
+                sprintf('Subject %d has no values in one or more conditions, setting all values to NaN', ...
+                nan_subs(z));
+            warning(msg);
+            project_amps(nan_subs(z),:,:) = NaN;
+        end
+    else
+    end
+     % make new valset and error set
+    valSet = squeeze(nanmean(project_amps,1));
+    temp_SE = squeeze(nanstd(project_amps,0,1)./sqrt(size(project_amps,1)));
+    errSet1 = temp_SE;
+    errSet2 = temp_SE;
+    % compute t-values
+    [~,temp_p,~,temp_stats] = ttest(project_amps,0,'alpha',0.05,'dim',1,'tail','right');
+    rc_tSqrdP(:,(1:2)+(f-1)*3) = permute(temp_p,[2,3,1]);
+    rc_tSqrdVal(:,(1:2)+(f-1)*3) = permute(temp_stats.tstat,[2,3,1]);
+    rc_tSqrdDF(:,(1:2)+(f-1)*3) = permute(temp_stats.df,[2,3,1]);
+    clear temp_*;
+
+    % maximum noise across all four conditions being plotted
+    % since this is just means, we can compute it the same way for
+    % projected and not-projected
+    noiseSet = max(superRCA(f).stats.NoiseAmp(:,rcNum,:),[],3);
+    % max of mean over bins
+    noiseSet(length(binVals)+1) = max(sqrt(nanmean(realBinMeanNoise,1).^2+nanmean(imagBinMeanNoise,1).^2));
+
+    for b = 1:(length(binVals)+1)
+        curData = squeeze(project_amps(:,b,:));
+        not_nan = ~any(isnan(curData),2);
+        [rc_tSqrdSig(b,3+(f-1)*3),rc_tSqrdP(b,3+(f-1)*3),~,tempStrct] = ttest(curData(not_nan,1),curData(not_nan,2),'alpha',0.05,'dim',1,'tail','both');
+        rc_tSqrdVal(b,3+(f-1)*3) = tempStrct.tstat;
+        rc_tSqrdDF(b,3+(f-1)*3) = tempStrct.df;
+        rc_tSqrdD(b,f) = rc_tSqrdVal(b,3+(f-1)*3)./sqrt(tempStrct.df+1); % Cohen's D
+        OVL = 2*normcdf(-abs(rc_tSqrdD(b,f)/2));
+        rc_tSqrdU(b,f) = 1-OVL/(2-OVL); 
+        rc_tSqrdMu1(b,f) = mean(curData(not_nan,1));
+        rc_tSqrdMu2(b,f) = mean(curData(not_nan,2));
+    end
+
+    % do Naka-Rushton, paired tests
+    testVal = squeeze(superRCA(f).stats.NR_JKParams(:,:,rcNum,:));
+    testVal = permute(testVal,[2,1,3]); % move subjects to first dim    
+    paramIdx = [1,2,3,4];% only look at c50 and rMax
+    jkDf = size(testVal,1)-1;
+    diffErr = jackKnifeErr(testVal(:,paramIdx,1)-testVal(:,paramIdx,2));
+    grandDiff = NRset(paramIdx,1) - NRset(paramIdx,2);
+    paramPairedT(:,f) = grandDiff'./diffErr;
+    paramPairedP(:,f) = 2*tcdf( -abs(paramPairedT(:,f)) , jkDf);
+
+    egiH(f) = subplot(2,2,2+(f-1)*2);
+    hold on
+    if f == 1
+        [figH(f),cH(f)] = mrC.plotOnEgi(superRCA(f).A(:,rcNum)*flipIdx(f),rcaColorBar,true);
+    else
+        [figH(f),cH(f)] = mrC.plotOnEgi(superRCA(f).A(:,rcNum)*flipIdx(f),rcaColorBar,true);
+    end
+    hold off
+    figH(f+2) = subplot(2,2,f+(f-1));
+    hold on
+    for c = 1:length(condsToUse)
+        if c == 1
+            ampMarkerStyle = {'o','LineWidth',lWidth,'Color',mainColors(c,:),'markerfacecolor',[1 1 1],'MarkerSize',7}; % 'MarkerEdgeColor','none'
+        else
+            ampMarkerStyle = {'sq','LineWidth',lWidth,'Color',mainColors(c,:),'markerfacecolor',[1 1 1],'MarkerSize',7}; % 'MarkerEdgeColor','none'
+        end
+        xVals = reallog([binVals;extraBins((condsToUse(c)>4)+1)]);
+        yVals = valSet(:,c);
+        errVals = [errSet1(:,c),errSet2(:,c)];
+        valH(c)=plot(xVals,yVals,ampMarkerStyle{:});
+        hE = ErrorBars(xVals,yVals,errVals,'color',mainColors(c,:),'type','bar','cap',false,'barwidth',lWidth);
+        uistack(valH(c),'bottom')
+        cellfun(@(x) uistack(x,'bottom'), hE);
+        hold on
+
+        % plot Naka-Rushton
+        nFine = 1e2;
+        nrX = linspace( min(binVals), max(binVals), nFine )';
+        nrVals = NRmodel( nrX, NRset(:,c));
+        hNR{c} = plot( reallog(nrX), nrVals, '-','color',mainColors(c,:), 'LineWidth',lWidth);
+    end
+    cellfun(@(x) uistack(x,'bottom'), hNR);
+
+    if f == 1     
+        yUnit = 1;
+        yMin = 0;
+        yMax = 5;
+        legend(valH,{'horizontal','vertical'},'fontsize',fSize,'fontname','Helvetica','location','northwest');
+        legend boxoff
+        titleStr = '\it\fontname{Helvetica}2F ';
+        title(titleStr,'fontsize',fSize,'fontname','Helvetica','interpreter','tex');
+    else
+        yUnit = 0.5;
+        yMin = 0;
+        yMax = 1.5;
+        ylabel('amplitude (\muV)','fontsize',fSize,'fontname','Helvetica')
+        xlabel('displacement (arcmins)','fontsize',fSize,'fontname','Helvetica');
+        titleStr = '\it\fontname{Arial}4F ';
+        title(titleStr,'fontsize',fSize,'fontname','Helvetica','interpreter','tex');
+    end
+    set(gca,gcaOpts{:},'XMinorTick','off','xtick',reallog([0.2,0.5,1,2,4,8,16]),'xticklabels',arrayfun(@(x) num2str(x),([0.2,0.5,1,2,4,8,16]),'uni',false),'ytick',0:yUnit:yMax,'Layer','top');
+
+    xlim([xMin,xMax]);
+    ylim([yMin,yMax])
+    text(xMin-logStep*.75,max(get(gca,'ylim'))+diff(get(gca,'ylim'))*.1,...
+        plotLabel(f),'fontsize',fSize*2,'fontname','Helvetica');
+    text(xMax+logStep*1.5,max(get(gca,'ylim'))+diff(get(gca,'ylim'))*.1,...
+        plotLabel(f+2),'fontsize',fSize*2,'fontname','Helvetica');
+
+    % split point between bin and average vals
+    xSplit = reallog(binVals(end))+logStep*.5;
+    plot(ones(2,1)*xSplit,[0,yMax],'k','LineWidth',lWidth)
+    % plot noise patch, 
+    % add mean bin noise twice to the end with extra bin
+    % we are plotting the mean values over two x-axis points
+    xNoiseVals = [xMin,xMin,reallog(binVals'),xSplit,xSplit];
+    yNoiseVals = [0,noiseSet(1),noiseSet(1:end-1)',noiseSet(end-1),0]; % start and end points just repeats of first and last
+    pH(1) = fill(xNoiseVals,yNoiseVals,[.75 .75 .75],'edgecolor','none');
+    % additional vals for averages
+    xNoiseVals = [xSplit,xSplit,reallog(extraBins),xMax,xMax];
+    yNoiseVals = [0,repmat(noiseSet(end),1,4),0];
+    pH(2) = fill(xNoiseVals,yNoiseVals,[.75 .75 .75],'edgecolor','none');
+    arrayfun(@(x) uistack(x,'bottom'), pH);
+end
+
+for f = 1:length(freqsToUse)
+    addX = 0.38;
+    addY = 0.38;
+    if f == length(freqsToUse)
+        set(cH(f),'fontsize',fSize,'fontname','Helvetica','YTick',linspace(min(rcaColorBar),min(rcaColorBar)*-1,5));
+        ylabel(cH(f),'weights','fontsize',fSize,'fontname','Helvetica')
+        set(cH(f),'location','eastoutside');
+        set(cH(f),'units','centimeters');
+        %cBarPos = get(cH(e,f),'position');
+        cBarPos = [14,4.5,.4,4];
+        set(cH(f),'position',cBarPos);
+        oldPos = newPos;
+        newPos = get(egiH(f),'position');
+        newPos(1) = oldPos(1);
+        newPos(2) = newPos(2)-(newPos(4)*addY*.6);
+        newPos(3) = oldPos(3);
+        newPos(4) = oldPos(4);
+    else
+        set(cH(f),'visible','off');
+        newPos = get(egiH(f),'position');
+        newPos(1) = newPos(1)-(newPos(3)*addX);
+        newPos(2) = newPos(2)-(newPos(4)*addY*.6);
+        newPos(3) = newPos(3)*(1+addX);
+        newPos(4) = newPos(4)*(1+addY);
+    end
+    set(egiH(f),'position',newPos); 
+end
+drawnow;
+set(gcf, 'units', 'centimeters');
+figPos = get(gcf,'pos');
+figPos(4) = figPos(4)/figPos(3)*17.8;
+figPos(3) = 17.8;
+set(gcf,'pos',figPos);
+
+