% 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