
% axesMin=[1, 20, 80, 0, 0, 0, 0,0,0, 250,50,50,220,80];
% axesMax=[3, 80, 100, 900, 900, 900, 2.5,2.5, 2.5, 330,120,120,320,100];
% axesLabels={'Group', 'Age (years)', 'Letter Score', 'V1 surface area (mm^2)', 'V2 surface area (mm^2)', 'V3 surface area (mm^2)', 'V1 global pRF size (?, °)', 'V2 global pRF size (?, °)', 'V3 global pRF size (?, °)',...
% 'Retinal thickness (?m)', 'Ganglion cell inner plexiform layer thickness (?m)', 'Retinal nerve fiber layer thickness (?m)', 'Central subfield thickness (?m)'};
% axesLabels={'Group', 'Age (years)', 'Letter Score', 'V1 surface area (mm2)', 'V2 surface area (mm2)', 'V3 surface area (mm2)', 'V1 global pRF size ()', 'V2 global pRF size ()', 'V3 global pRF size ()',...
% 'Retinal thickness (um)', 'Ganglion cell inner plexiform layer thickness (um)', 'Retinal nerve fiber layer thickness (um)', 'Central subfield thickness (m)'};
% 
% %Some are measures of eyes, some are hemispheres. These can't be mixed, so
% %are averages across eyes where mixed.
% eyeList=[2 10 11 12 13];
% hemiList=[2 4 5 6 7 8 9];
% 
% %To rearrange orders (optional)
% %newOrder=[1 3 7 8 9 4 5 6 10 11 12 13, 2];%letter score, pRFs, Areas, RT, age
newOrder=[10 11 12 4 5 6 7 8 9 3 2]; %For correlation grids
AgingMeans=AgingMeans(:,newOrder); 
Aging2Eyes=Aging2Eyes(:,newOrder);
% %axesLabels=axesLabels{newOrder};
% for n=1:length(newOrder)
%     axesLabelsNew{n}=axesLabels{newOrder(n)};
% end
% axesLabels=axesLabelsNew;
% [~,eyeList]=intersect(newOrder,eyeList);
% [~,hemiList]=intersect(newOrder,hemiList);
% axesMin=axesMin(newOrder);
% axesMax=axesMax(newOrder);

axesMin=[250 50 50 0 0 0 0 0 0 80 20];
axesMax=[330 120 120 900 900 900 2.5 2.5 2.5 100 80];
eyeList=[1 2 3 11]; %Measures per eye
hemiList=[4 5 6 7 8 9 11]; %Measures per hemisphere
axesLabels={'Retinal thickness (um)', 'Ganglion cell inner plexiform layer thickness (um)', 'Retinal nerve fiber layer thickness (um)', 'V1 surface area (mm2)', 'V2 surface area (mm2)', 'V3 surface area (mm2)', 'V1 global pRF size ()', 'V2 global pRF size ()', 'V3 global pRF size ()','Letter Score','Age (years)'}; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Only run the part above once.
%Choose which parameters to correlate. 
%1= Retinal thickness,2=GCIPL thickness, 3=RNFL thickness, 4,5,6=Surface area V1 V2 V3; 7,8,9=pRF size V1 V2 V3; 10=Corrected visual acuity (letter score) 11=Age
Param1=4;
Param2=1;

%For one specific correlation plot
if (ismember(Param1, eyeList) && ismember(Param2, eyeList)) || (ismember(Param1, hemiList) && ismember(Param2, hemiList))
    [r,p]=corr(Aging2Eyes(~isnan(Aging2Eyes(:,Param1)+Aging2Eyes(:,Param2)),Param1), Aging2Eyes(~isnan(Aging2Eyes(:,Param1)+Aging2Eyes(:,Param2)),Param2), 'Type', 'Spearman');
    [fit] = polyfit(Aging2Eyes(~isnan(Aging2Eyes(:,Param1)+Aging2Eyes(:,Param2)),Param1), Aging2Eyes(~isnan(Aging2Eyes(:,Param1)+Aging2Eyes(:,Param2)),Param2),1);
    xs=[min(Aging2Eyes(:,Param1)) max(Aging2Eyes(:,Param1))];
    ys=polyval(fit, xs);
    figure; plot(Aging2Eyes(Aging2Eyes(:,11)<40,Param1), Aging2Eyes(Aging2Eyes(:,11)<40,Param2), 'ro');
    hold on; plot(Aging2Eyes((Aging2Eyes(:,11)>=40 & Aging2Eyes(:,11)<60),Param1), Aging2Eyes((Aging2Eyes(:,11)>=40 & Aging2Eyes(:,11)<60),Param2), 'go');
    hold on; plot(Aging2Eyes((Aging2Eyes(:,11)>=60 & Aging2Eyes(:,11)<80),Param1), Aging2Eyes((Aging2Eyes(:,11)>=60 & Aging2Eyes(:,11)<80),Param2), 'bo');
    hold on; plot(xs,ys, 'k')
    n=100;
else %If one is an eye measure and the other is hemisphere, use the average.
    [r,p]=corr(AgingMeans(~isnan(AgingMeans(:,Param1)+AgingMeans(:,Param2)),Param1), AgingMeans(~isnan(AgingMeans(:,Param1)+AgingMeans(:,Param2)),Param2), 'Type', 'Spearman');
    [fit] = polyfit(AgingMeans(~isnan(AgingMeans(:,Param1)+AgingMeans(:,Param2)),Param1), AgingMeans(~isnan(AgingMeans(:,Param1)+AgingMeans(:,Param2)),Param2),1);
    xs=[min(AgingMeans(:,Param1)) max(AgingMeans(:,Param1))];
    ys=polyval(fit, xs);
    figure; plot(AgingMeans(AgingMeans(:,11)<40,Param1), AgingMeans(AgingMeans(:,11)<40,Param2), 'ro');
    hold on; plot(AgingMeans((AgingMeans(:,11)>=40 & AgingMeans(:,11)<60),Param1), AgingMeans((AgingMeans(:,11)>=40 & AgingMeans(:,11)<60),Param2), 'go');
    hold on; plot(AgingMeans((AgingMeans(:,11)>=60 & AgingMeans(:,11)<80),Param1), AgingMeans((AgingMeans(:,11)>=60 & AgingMeans(:,11)<80),Param2), 'bo');
    hold on; plot(xs,ys, 'k')
    n=50;
end
axis square;
xlabel(axesLabels{Param1});
ylabel(axesLabels{Param2});
title(sprintf('r = %.2f, n = %.0f, p = %.4f', r,n,p))
fit(1)
axis([axesMin(Param1), axesMax(Param1), axesMin(Param2), axesMax(Param2)])




%%%%%%%%%%%
%To make a grid showing all correlations.
pall=[];
nall=[];
rall=[];

for Param1=1:11
    for Param2=1:11
        if (ismember(Param1, eyeList) && ismember(Param2, eyeList)) || (ismember(Param1, hemiList) && ismember(Param2, hemiList))
            [rall(Param1, Param2), pall(Param1, Param2)]=corr(Aging2Eyes(~isnan(Aging2Eyes(:,Param1)+Aging2Eyes(:,Param2)),Param1), Aging2Eyes(~isnan(Aging2Eyes(:,Param1)+Aging2Eyes(:,Param2)),Param2), 'Type', 'Spearman');
            nall(Param1, Param2)=100;
        else
            [rall(Param1, Param2), pall(Param1, Param2)]=corr(AgingMeans(~isnan(AgingMeans(:,Param1)+AgingMeans(:,Param2)),Param1), AgingMeans(~isnan(AgingMeans(:,Param1)+AgingMeans(:,Param2)),Param2), 'Type', 'Spearman');
            nall(Param1, Param2)=50;
       end
    end
end

%False discovery rate correction for multiple comparisons.
[~,~,~,pall]=fdr_bh(pall);

%Plot correlation grid
logP=-log10(pall);
logP(rall<0)=-logP(rall<0);

rallBig=[];
logPBig=[];
for Param1=1:11
    for Param2=1:11
        rallBig(((Param1-1)*10+1):(Param1*10), ((Param2-1)*10+1):(Param2*10))=rall(Param1,Param2);
        logPBig(((Param1-1)*10+1):(Param1*10), ((Param2-1)*10+1):(Param2*10))=logP(Param1,Param2);
    end
end
figure; imagesc(rallBig, [-1 1])
axis image;
colorbar
colormap jet
for Param1=1:11
        text((Param1-1)*10+5,1,axesLabels{Param1})
        text(1,(Param1-1)*10+5,axesLabels{Param1})
end

figure; imagesc(logPBig, [-4 4])
axis image;
colorbar
colormap jet
for Param1=1:11
        text((Param1-1)*10+5,1,axesLabels{Param1})
        text(1,(Param1-1)*10+5,axesLabels{Param1})
end

figure; imagesc(abs(logPBig), [0 6])
axis image;
colorbar
colormap hot
for Param1=1:11
        text((Param1-1)*10+5,1,axesLabels{Param1})
        text(1,(Param1-1)*10+5,axesLabels{Param1})
end


%Bayes factor image
%Taken from JASP
logBF=[0 18.963 6.752 1.549 3.512 2.076 -0.516 -1.344 -1.446 -0.866 1.162;...
    18.963 0 4.311 2.875 4.441 3.656 0.365 -.443 -0.738 1.491 0.880;...
    6.752 4.311 0 2.129 2.361 1.149 -1.265 -1.527 -0.351 0.631 1.130;...
    1.549 2.875 2.129 0 15.251 13.120 -0.340 -0.830 -0.172 -0.051 2.707;...
    3.512 4.441 2.361 15.251 0 18.344 0.743 1.044 -0.599 -1.134 3.754;...
    2.076 3.656 1.149 13.120 18.344 0 0.099 0.292 -0.143 -1.672 1.555;...
    -0.516 0.365 -1.265 -0.340 0.743 0.099 0 26.614 13.390 2.960 5.157;...
    -1.344 -0.443 -1.527 -0.830 1.044 0.292 26.614 0 15.251 1.285 3.217;...
    -1.446 -0.738 -0.351 -0.172 -0.599 -0.143 13.390 15.251 0 1.930 4.494;...
    -0.866 -1.491 0.631 -0.051 -1.134 -1.672 2.960 1.285 1.930 0 10.866;...]
    1.162 0.880 1.130 2.707 3.754 1.555 5.157 3.217 4.494 10.866 0];

logBF100=[0	40.732	12.891	5.453	4.502	2.069	1.617	-1.279	-1.242	-0.11	4.633;...
40.732	0	10.135	5.498	6.147	5.182	1.82	0.425	-0.4	-1.352	3.736;...
12.891	10.135	0	0.885	2.312	1.002	-1.288	-1.918	-0.47	1.88	3.29;...
5.453	5.498	0.885	0	26.204	14.876	0.927	-1.1	-0.661	1.204	5.777;...
4.502	6.147	2.312	26.204	0	32.358	3.127	3.127	-1.508	-1.111	5.184;...
2.069	5.182	1.002	14.876	32.358	0	2.386	1.645	-0.611	-1.966	2.528;...
1.617	1.82	-1.288	0.927	3.127	2.386	0	35.849	16.85	7.026	10.824;...
-1.279	0.425	-1.918	-1.1	3.127	1.645	35.849	0	26.31	3.119	7.062;...
-1.242	-0.4	-0.47	-0.661	-1.508	-0.611	16.85	26.31	0	3.456	7.451;...
-0.11	-1.352	1.88	1.204	-1.111	-1.966	7.026	3.119	3.456	0	24.286;...
4.633	3.736	3.29	5.777	5.184	2.528	10.824	7.062	7.451	24.286	0];

logBF(nall==100)=logBF100(nall==100);
logBF(logBF==0)=5;

%Convert to log10
BF=exp(logBF);
logBF=log10(BF);

logBFBig=[];
for Param1=1:11
    for Param2=1:11
        logBFBig(((Param1-1)*10+1):(Param1*10), ((Param2-1)*10+1):(Param2*10))=logBF(Param1,Param2);
    end
end
figure; imagesc(logBFBig, [-2 2])
axis image;
cbh=colorbar;
cbh.Ticks=log10([0.01 0.02 1/20 0.1 0.2 1/2 1 2 5 10 20 50 100])
colormap jet
for Param1=1:11
        text((Param1-1)*10+5,1,axesLabels{Param1})
        text(1,(Param1-1)*10+5,axesLabels{Param1})
end


%%%%%%%%%%%%
%Mediation models
%1 is the effect of age on acuity mediated by V1 pRF size?
Param1=12;
p1=axesLabels{Param1}
Param2=5;
p2=axesLabels{Param2}
Mediator=2;
med=axesLabels{Mediator}
X=AgingMeans(:,Param1);
Y=AgingMeans(:,Param2);
M=AgingMeans(:,Mediator);
mask=~isnan(AgingMeans(:,Param1)+AgingMeans(:,Param2)+AgingMeans(:,Mediator));

% 1 a   X -> M relationship
% 2 b   M -> Y relationship
% 3 cp  unmediated X -> Y relationship (residual)
% 4 c   X -> Y relationship
% 5 ab  mediated X -> Y by M (a * b)

[paths, toplevelstats] = mediation(X, Y, M, 'names', {p1 p2 med})


%GLM
allData(:,1)=AgingMeans(:,1); %RT
allData(:,2)=AgingMeans(:,4); %Area 
allData(:,3)=AgingMeans(:,7); %pRF
allData(:,4)=AgingMeans(:,10); %Acuity
age=AgingMeans(:,11);

allData=allData([1:28 30:end],:);%Exclude missing data

allT=nan(4,4);
for x=1:4
    whichX=1:4;
    whichX=whichX~=x;
    X=[allData(:,whichX) ones(size(allData(:,1)))];
    Y=allData(:,x);
    B_hat=pinv(X)*Y;
    U=Y-X*B_hat;
    df = size(Y,1) - size(X,2);
    for y=1:3
            C=[0 0 0 0];
            C(y)=1;
            SE = sqrt((sum(U.^2)./df)*(C*pinv(X'*X)*C')); 
            if y<x
                allT(x,y)=C*B_hat./SE;
            else
               allT(x,y+1)=C*B_hat./SE; 
            end
    end
end
allT(isnan(allT))=5;

allP=2*(1 - tpvalue(abs(allT), df));

rall=allT;
logP=-log10(allP);

rallBig=[];
logPBig=[];
for Param1=1:4
    for Param2=1:4
        rallBig(((Param1-1)*10+1):(Param1*10), ((Param2-1)*10+1):(Param2*10))=rall(Param1,Param2);
        logPBig(((Param1-1)*10+1):(Param1*10), ((Param2-1)*10+1):(Param2*10))=logP(Param1,Param2);
    end
end
figure; imagesc(rallBig, [-4 4])
axis image;
colorbar
colormap jet
whichLabels=[9 6 3 2];
for Param1=1:4
        text((Param1-1)*10+5,1,axesLabels{whichLabels(Param1)})
        text(1,(Param1-1)*10+5,axesLabels{whichLabels(Param1)})
end

figure; imagesc(logPBig, [0 4])
axis image;
colorbar
colormap hot
for Param1=1:4
        text((Param1-1)*10+5,1,axesLabels{whichLabels(Param1)})
        text(1,(Param1-1)*10+5,axesLabels{whichLabels(Param1)})
end



