From b2e1f9d4e0a868157dd8f76de9de45cd741c113a Mon Sep 17 00:00:00 2001 From: Atsushi Yokoi Date: Tue, 6 Oct 2015 19:56:21 -0400 Subject: [PATCH 1/3] Updated sh1_rsa - added a case for group reliability - added a case for estimating tau with pooling ROI data --- sh1_rsa.m | 338 +++++++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 312 insertions(+), 26 deletions(-) diff --git a/sh1_rsa.m b/sh1_rsa.m index 06cb00e..352c17a 100644 --- a/sh1_rsa.m +++ b/sh1_rsa.m @@ -229,6 +229,65 @@ end; varargout={T}; save(fullfile(regDir,'distances_sepPerm.mat'),'-struct','T'); + case 'ROI_reliability_btw' % Calculate between subject reliability of RDMs + + C = []; + for glm = [1 3] + % load data set + T = load(fullfile(regDir,sprintf('distances_sepPerm.glm%d.mat',glm))); + + chunksets = {'A','B'}; + for reg = [1:8] + for h = [1 2] + for chunkset = [1 2] + % calc within chunkset value + D = getrow(T,T.region==reg&T.hemis==h&T.chunkset==chunksets{chunkset}); + D_ = getrow(T,T.region==reg&T.hemis==h&T.chunkset~=chunksets{chunkset}); + + rdm_all = nanmean(D.RDM,1); % all rdm for the same chunkset + rdm_all_ = nanmean(D_.RDM,1); % all rdm for other chunkset + subj = unique(D.subj); + for s=1:numel(subj) + idx = D.subj==subj(s); + rdm = D.RDM(idx,:); + rdm_loo = nanmean(D.RDM(~idx,:),1); % leave-one-subject-out rdm + + R.within_u = corr(rdm',rdm_all','type','Pearson'); % upper bound + R.within_l = corr(rdm',rdm_loo','type','Pearson'); % lower bound + + % we must choose only the distance between physically identical sequences here + tmprdm = getCommon(rdm,[3 5 7]); + tmprdm_ = getCommon(rdm_all_,[3 5 7]); + R.across = corr(tmprdm',tmprdm_','type','Pearson'); % across chunk set corr + + R.subj = subj(s); + R.chunkset = chunkset; + R.region = reg; + R.hemis = h; + R.glm = glm; + + C = addstruct(C,R); + end + end + end + end + end + + for glm=[1 3] + figure('name',sprintf('reliability (glm=%d)',glm)); + for reg=1:8; + subplot(2,4,reg); + %xpos = myboxplot(C.hemis,[C.within_u,C.within_l,C.across],'subset',C.region==reg); + xpos = barplot([C.hemis],[C.within_u,C.within_l,C.across],'subset',C.region==reg&C.glm==glm); + drawline(0,'dir','horz'); + title(regname{reg},'fontsize',12); + ylabel('Pearson r','fontsize',12) + %set(gca,'XTick',[xpos(1:4)],'XTicklabel',omega_names,'fontsize',12); + %rotateXLabels(gca,45); + end + end + + varargout = {C}; case 'fit_model_EB_lin' T=load(fullfile(regDir,'distances_sepPerm.mat')); @@ -293,7 +352,9 @@ Model.nonlinP0 = [-10 1 10 0];%zeros(size(Model.constantParams{2})); % Starting value of nonlinear(mixing) parameters % estimate - [w,logEvidence(indx,:),lt]=rsa_fitModelHierarchEB(Model,D.RDM,Sigma,9,D.effVox); + [w,logEvidence(indx,:),lt] = rsa_fitModelHierarchEB(Model,D.RDM,Sigma,9,D.effVox); + + % logtheta(indx,:)=repmat(lt,length(indx),1); % normalize omega with mean of each RDM @@ -321,14 +382,165 @@ T.logEvidence=logEvidence; %T.logEvidenceSplit = logEvidenceSplit; - sh1_rsa('plot_omega',T); + sh1_rsa('plot_omega_prior',T); + sh1_rsa('plot_prior',T); + sh1_rsa('plot_tau',T); + varargout={T,mRDM}; + case 'fit_model_EB_nonlin_poolROI' + glm = 3; + useolddata = 0; + modelTerms = [1 2 4 6]; + + vararginoptions(varargin(:),{'glm','useolddata','modelTerms'}) + + % load data set + if useolddata + %T = load(fullfile(regDir,'distances_glm1.mat')); + T = load(fullfile(regDir,'allRDMs_ROI_all_glm3.mat')); + %T = T.RDM1; + T.RDM = bsxfun(@rdivide,T.RDM,T.numVox); + else + T = load(fullfile(regDir,sprintf('distances_sepPerm.glm%d.mat',glm))); + end + + % pool all ROI data to estimate single tau parameter for whole brain + mRDM = []; + indx = true(size(T.subj));%find(T.region==1&T.hemis==1); + D=getrow(T,indx); + for s=1:length(D.subj) + Sigma(:,:,s)=reshape(D.Sigma(s,:),8,8); + end; + + for iter = 1:5 + iter + initialguess = rand(size(modelTerms))*100-50; + initialguess(end) = 0; + + % organize model RDMs + Model.fcn = @sh1_getRDMmodelTau1; + Model.constantParams = {chunk_set(D.subj),... % Model terms + modelTerms,... % Chunk set + 'sqEuclidean'... % Distance term + }; + Model.numComp = numel(Model.constantParams{2}); % Number of linear components + Model.numPrior = numel(Model.constantParams{2}); % Number of prior variances on parameters + Model.numNonlin = numel(Model.constantParams{2}); % Number of nonlinear parameters + Model.nonlinP0 = initialguess; % Starting value of nonlinear(mixing) parameters + + % estimate + [w(:,:,iter),levid(:,iter),lt(iter,:)] = rsa_fitModelHierarchEB(Model,D.RDM,Sigma,9,D.effVox); + + mevid = mean(levid,1); + end + maxidx = find(mevid==max(mevid),1,'first'); + % choose one + logtheta(indx,:) = repmat(lt(maxidx,:),length(indx),1); + w = squeeze(w(:,:,maxidx)); + logEvidence(indx,:) = levid(:,maxidx); + + % normalize omega with mean of each RDM + M = feval(Model.fcn,lt(Model.numPrior+1:Model.numPrior+Model.numNonlin),... + Model.constantParams{:}); + norm = permute(mean(M.RDM.^2,2),[3 1 2]); + omega(indx,:) = w; + omegaN(indx,:) = w.*norm; + + % save resultant RDMs with estimated tau +% A = feval(Model.fcn,lt(Model.numPrior+1:Model.numPrior+Model.numNonlin),... +% 1,Model.constantParams{2:end}); % chunk set A +% B = feval(Model.fcn,lt(Model.numPrior+1:Model.numPrior+Model.numNonlin),... +% 2,Model.constantParams{2:end}); % chunk set B +% mrdm.RDM = [A.RDM;B.RDM]; +% mrdm.set = kron([1 2]',ones(Model.numComp,1)); +% mrdm.comp = repmat(Model.constantParams{2}',2,1); +% mrdm.name = {A.name{:};B.name{:}}; +% mrdm.region = repmat(r,size(mrdm.set)); +% mrdm.hemis = repmat(h,size(mrdm.set)); +% mRDM = addstruct(mRDM,mrdm); +% % + + T.omega = omega; + T.omegaN = omegaN; + T.logtheta = logtheta; + T.logEvidence = logEvidence; + %T.logEvidenceSplit = logEvidenceSplit; + + sh1_rsa('plot_omega_prior',T); sh1_rsa('plot_prior',T); sh1_rsa('plot_tau',T); varargout={T,mRDM}; case 'fit_model_EB_nonlin_singlemodel' % fit single model to compare log-evidence + glm = 1; + useolddata = 0; + modelTerms = [1 2 4 6]; + nonlinP0 = [-10 1 10 0]; + vararginoptions(varargin(:),{'glm','useolddata','modelTerms','nonlinP0'}) + + % load data set + if useolddata + %T = load(fullfile(regDir,'distances_glm1.mat')); + T = load(fullfile(regDir,'allRDMs_ROI_all_glm3.mat')); + %T = T.RDM1; + T.RDM = bsxfun(@rdivide,T.RDM,T.numVox); + else + T = load(fullfile(regDir,sprintf('distances_sepPerm.glm%d.mat',glm))); + end + + %regions=unique(T.region); + regions=[1:8]'; + for r=regions' + fprintf('Estimating ROI distance %s...\n',regname{r}); + for h=[1 2] + for m=1:numel(modelTerms) + indx = find(T.region==r & T.hemis==h); + D=getrow(T,indx); + for s=1:length(D.subj) + Sigma(:,:,s)=reshape(D.Sigma(s,:),8,8); + end; + + % organize model RDMs + Model.fcn = @sh1_getRDMmodelTau1; + Model.constantParams = {chunk_set(D.subj),... % Model terms + modelTerms(m),... % Chunk set + 'sqEuclidean'... % Distance term + }; + Model.numComp = numel(Model.constantParams{2}); % Number of linear components + Model.numPrior = numel(Model.constantParams{2}); % Number of prior variances on parameters + Model.numNonlin = numel(Model.constantParams{2}); % Number of nonlinear parameters + Model.nonlinP0 = nonlinP0(m);%zeros(size(Model.constantParams{2})); % Starting value of nonlinear(mixing) parameters + + % estimate + [w,logEvidence(indx,m),lt] = rsa_fitModelHierarchEB(Model,D.RDM,Sigma,9,D.effVox); + + logprior(indx,m) = repmat(lt(1),length(indx),1); + logtau(indx,m) = repmat(lt(2),length(indx),1); + term(indx,m) = repmat(modelTerms(m),length(indx),1);; + + % normalize omega with mean of each RDM + M = feval(Model.fcn,lt(Model.numPrior+1:Model.numPrior+Model.numNonlin),... + Model.constantParams{:}); + omega(indx,m) = w.*permute(mean(M.RDM,2),[3 1 2]); + + end % model term + end; + end; + T.omega = omega; + T.logtheta = [logprior, logtau]; + T.logEvidence = logEvidence; + T.term = term; + %T.logEvidenceSplit = logEvidenceSplit; + + % plot result + sh1_rsa('plot_logEvidence',T); + sh1_rsa('plot_omega',T); + sh1_rsa('plot_prior',T); + sh1_rsa('plot_tau',T); + varargout={T}; + case 'fit_model_EB_nonlin_multiguess' glm = 1; useolddata = 1; - vararginoptions(varargin(:),{'glm','useolddata'}) + modelTerms = [1 2 4 6]; + vararginoptions(varargin(:),{'glm','useolddata','modelTerms'}) % load data set if useolddata @@ -341,7 +553,8 @@ end %regions=unique(T.region); - regions=[1:8]'; + regions=[1:8]'; + mRDM = []; for r=regions' for h=[1 2] indx = find(T.region==r & T.hemis==h); @@ -350,29 +563,51 @@ Sigma(:,:,s)=reshape(D.Sigma(s,:),8,8); end; - % for model term - - % organize model RDMs - Model.fcn = @sh1_getRDMmodelTau1; - Model.constantParams = {chunk_set(D.subj),... % Model terms - [1 2 4 6],... % Chunk set - 'sqEuclidean'... % Distance term - }; - Model.numComp = numel(Model.constantParams{2}); % Number of linear components - Model.numPrior = numel(Model.constantParams{2}); % Number of prior variances on parameters - Model.numNonlin = numel(Model.constantParams{2}); % Number of nonlinear parameters - Model.nonlinP0 = [-10 1 10 0];%zeros(size(Model.constantParams{2})); % Starting value of nonlinear(mixing) parameters - - % estimate - [w,logEvidence(indx,:),lt]=rsa_fitModelHierarchEB(Model,D.RDM,Sigma,9,D.effVox); - logtheta(indx,:)=repmat(lt,length(indx),1); + for iter = 1:50 + iter + initialguess = rand(size(modelTerms))*100-50; + initialguess(end) = 0; + + % organize model RDMs + Model.fcn = @sh1_getRDMmodelTau1; + Model.constantParams = {chunk_set(D.subj),... % Model terms + modelTerms,... % Chunk set + 'sqEuclidean'... % Distance term + }; + Model.numComp = numel(Model.constantParams{2}); % Number of linear components + Model.numPrior = numel(Model.constantParams{2}); % Number of prior variances on parameters + Model.numNonlin = numel(Model.constantParams{2}); % Number of nonlinear parameters + Model.nonlinP0 = initialguess; % Starting value of nonlinear(mixing) parameters + + % estimate + [w(:,:,iter),levid(:,iter),lt(iter,:)] = rsa_fitModelHierarchEB(Model,D.RDM,Sigma,9,D.effVox); + + mevid = mean(levid,1); + end + maxidx = find(mevid==max(mevid),1,'first'); + % choose one + logtheta(indx,:) = repmat(lt(maxidx,:),length(indx),1); + w = squeeze(w(:,:,maxidx)); + logEvidence(indx,:) = levid(:,maxidx); % normalize omega with mean of each RDM M = feval(Model.fcn,lt(Model.numPrior+1:Model.numPrior+Model.numNonlin),... Model.constantParams{:}); - omega(indx,:) = w.*permute(mean(M.RDM,2),[3 1 2]); + omega(indx,:) = w.*permute(mean(M.RDM,2),[3 1 2]); - % end model term + % save resultant RDMs with estimated tau + A = feval(Model.fcn,lt(Model.numPrior+1:Model.numPrior+Model.numNonlin),... + 1,Model.constantParams{2:end}); % chunk set A + B = feval(Model.fcn,lt(Model.numPrior+1:Model.numPrior+Model.numNonlin),... + 2,Model.constantParams{2:end}); % chunk set B + mrdm.RDM = [A.RDM;B.RDM]; + mrdm.set = kron([1 2]',ones(Model.numComp,1)); + mrdm.comp = repmat(Model.constantParams{2}',2,1); + mrdm.name = {A.name{:};B.name{:}}; + mrdm.region = repmat(r,size(mrdm.set)); + mrdm.hemis = repmat(h,size(mrdm.set)); + mRDM = addstruct(mRDM,mrdm); + % end; end; T.omega=omega; @@ -380,11 +615,11 @@ T.logEvidence=logEvidence; %T.logEvidenceSplit = logEvidenceSplit; - sh1_rsa('plot_omega',T); + sh1_rsa('plot_omega_prior',T); sh1_rsa('plot_prior',T); sh1_rsa('plot_tau',T); - varargout={T}; - + varargout={T,mRDM}; + case 'plot_omega' T=varargin{1}; @@ -423,7 +658,43 @@ title(regname{reg},'fontsize',12); %set(gca,'XTick',[xpos(1:4)],'XTicklabel',omega_names,'fontsize',12); %rotateXLabels(gca,45); - end + end + case 'plot_logEvidence' + T=varargin{1}; + + figure('name','logEvidence'); + for reg=1:8; + subplot(2,4,reg); + xpos=myboxplot(T.hemis,T.logEvidence,'subset',T.region==reg,... + 'CAT',CAT); + drawline(0,'dir','horz'); + title(regname{reg},'fontsize',12); + %set(gca,'XTick',[xpos(1:4)],'XTicklabel',omega_names,'fontsize',12); + %rotateXLabels(gca,45); + end + case 'plot_omega_prior' + T=varargin{1}; + + figure('name','omega'); + for reg=1:8; + subplot(2,8,reg); + xpos=myboxplot(T.hemis,T.omega,'subset',T.region==reg,... + 'CAT',CAT); + drawline(0,'dir','horz'); + title(regname{reg},'fontsize',12); + %set(gca,'XTick',[xpos(1:4)],'XTicklabel',omega_names,'fontsize',12); + %rotateXLabels(gca,45); + + % plot log prior + subplot(2,8,reg+8); + data = getrow(T,T.region==reg); + scale = max(max(abs(data.logtheta(:,1:4)))); + Ldata = getrow(T,T.region==reg&T.hemis==1); + Rdata = getrow(T,T.region==reg&T.hemis==2); + imagesc_rectangle([mean(Ldata.logtheta(:,1:4),1),mean(Rdata.logtheta(:,1:4),1)],... + 'scale',[-scale, scale]); + end + case 'showRDMs' T=varargin{1}; glm=varargin{2}; @@ -437,7 +708,22 @@ %saveas(h,fullfile(figDir,figName),'fig'); end end + otherwise disp('there is no such case.') end; + function RDMout = getCommon(RDMin, seqnum); + % RDMin : vectorised RDM + % RDMout: vectorised RDM + % seqnum: seq of no interest + + I = true(nSequence); + I(diag(diag(I))) = false; + I(seqnum,:) = false; + I(:,seqnum) = false; + + RDMout = RDMin(:,squareform(I)); + end + +end From 261f18c66f5b8365693a4e8324d2ab637cc20898 Mon Sep 17 00:00:00 2001 From: Atsushi Yokoi Date: Wed, 7 Oct 2015 16:22:57 -0400 Subject: [PATCH 2/3] added between-chunkset, between-subject reliability of RDM --- sh1_rsa.m | 183 ++++++++++++++++++++++++++++++++++++------------------ 1 file changed, 124 insertions(+), 59 deletions(-) diff --git a/sh1_rsa.m b/sh1_rsa.m index 352c17a..a8445b4 100644 --- a/sh1_rsa.m +++ b/sh1_rsa.m @@ -257,9 +257,13 @@ % we must choose only the distance between physically identical sequences here tmprdm = getCommon(rdm,[3 5 7]); - tmprdm_ = getCommon(rdm_all_,[3 5 7]); + tmprdm_ = getCommon(rdm_all_,[3 5 7]); R.across = corr(tmprdm',tmprdm_','type','Pearson'); % across chunk set corr + tmprdm = getCommon(rdm,[3 5 7]); + tmprdm_ = getCommon(rdm_loo,[3 5 7]); + R.within_l_reduced = corr(tmprdm',tmprdm_','type','Pearson'); % lower bound (calculated with reduced RDM) + R.subj = subj(s); R.chunkset = chunkset; R.region = reg; @@ -275,16 +279,8 @@ for glm=[1 3] figure('name',sprintf('reliability (glm=%d)',glm)); - for reg=1:8; - subplot(2,4,reg); - %xpos = myboxplot(C.hemis,[C.within_u,C.within_l,C.across],'subset',C.region==reg); - xpos = barplot([C.hemis],[C.within_u,C.within_l,C.across],'subset',C.region==reg&C.glm==glm); - drawline(0,'dir','horz'); - title(regname{reg},'fontsize',12); - ylabel('Pearson r','fontsize',12) - %set(gca,'XTick',[xpos(1:4)],'XTicklabel',omega_names,'fontsize',12); - %rotateXLabels(gca,45); - end + D = getrow(C,C.glm==glm); + sh1_rsa('plot_reliability_btw',D); end varargout = {C}; @@ -390,9 +386,12 @@ glm = 3; useolddata = 0; modelTerms = [1 2 4 6]; - - vararginoptions(varargin(:),{'glm','useolddata','modelTerms'}) + Niter = 100; % iteration from different initial parameter values + Ntake = 100; % how many best models you choose + region = [1 2]; + vararginoptions(varargin(:),{'glm','useolddata','modelTerms','Niter','Ntake','region'}) + TT=[]; % load data set if useolddata %T = load(fullfile(regDir,'distances_glm1.mat')); @@ -405,16 +404,16 @@ % pool all ROI data to estimate single tau parameter for whole brain mRDM = []; - indx = true(size(T.subj));%find(T.region==1&T.hemis==1); + indx = find(ismember(T.region,region)&T.hemis==1);%true(size(T.subj));%find(T.region==1&T.hemis==1); D=getrow(T,indx); for s=1:length(D.subj) Sigma(:,:,s)=reshape(D.Sigma(s,:),8,8); end; - for iter = 1:5 + for iter = 1:Niter iter - initialguess = rand(size(modelTerms))*100-50; - initialguess(end) = 0; + initialguess = rand(size(modelTerms))*3-1.5; + initialguess(end) = 0 % organize model RDMs Model.fcn = @sh1_getRDMmodelTau1; @@ -432,43 +431,65 @@ mevid = mean(levid,1); end - maxidx = find(mevid==max(mevid),1,'first'); - % choose one - logtheta(indx,:) = repmat(lt(maxidx,:),length(indx),1); - w = squeeze(w(:,:,maxidx)); - logEvidence(indx,:) = levid(:,maxidx); - - % normalize omega with mean of each RDM - M = feval(Model.fcn,lt(Model.numPrior+1:Model.numPrior+Model.numNonlin),... - Model.constantParams{:}); - norm = permute(mean(M.RDM.^2,2),[3 1 2]); - omega(indx,:) = w; - omegaN(indx,:) = w.*norm; - - % save resultant RDMs with estimated tau -% A = feval(Model.fcn,lt(Model.numPrior+1:Model.numPrior+Model.numNonlin),... -% 1,Model.constantParams{2:end}); % chunk set A -% B = feval(Model.fcn,lt(Model.numPrior+1:Model.numPrior+Model.numNonlin),... -% 2,Model.constantParams{2:end}); % chunk set B -% mrdm.RDM = [A.RDM;B.RDM]; -% mrdm.set = kron([1 2]',ones(Model.numComp,1)); -% mrdm.comp = repmat(Model.constantParams{2}',2,1); -% mrdm.name = {A.name{:};B.name{:}}; -% mrdm.region = repmat(r,size(mrdm.set)); -% mrdm.hemis = repmat(h,size(mrdm.set)); -% mRDM = addstruct(mRDM,mrdm); -% % - - T.omega = omega; - T.omegaN = omegaN; - T.logtheta = logtheta; - T.logEvidence = logEvidence; - %T.logEvidenceSplit = logEvidenceSplit; + [~,indices]= sort(mevid,2,'descend'); + maxindices = indices(1:Ntake);%find(mevid==max(mevid),Ntake,'first'); - sh1_rsa('plot_omega_prior',T); - sh1_rsa('plot_prior',T); - sh1_rsa('plot_tau',T); - varargout={T,mRDM}; + for m=1:length(maxindices) + + maxidx=maxindices(m); + % choose + logtheta(:,:) = repmat(lt(maxidx,:),length(indx),1); + ww = squeeze(w(:,:,maxidx)); + logEvidence(:,:) = levid(:,maxidx); + + % normalize omega with mean of each RDM + M = feval(Model.fcn,lt(Model.numPrior+1:Model.numPrior+Model.numNonlin),... + Model.constantParams{:}); + normX = sqrt(mean(M.RDM.^2,2)); + omega(:,:) = ww; + omegaN(:,:) = bsxfun(@times,ww,normX'); + + % save resultant RDMs with estimated tau + % A = feval(Model.fcn,lt(Model.numPrior+1:Model.numPrior+Model.numNonlin),... + % 1,Model.constantParams{2:end}); % chunk set A + % B = feval(Model.fcn,lt(Model.numPrior+1:Model.numPrior+Model.numNonlin),... + % 2,Model.constantParams{2:end}); % chunk set B + % mrdm.RDM = [A.RDM;B.RDM]; + % mrdm.set = kron([1 2]',ones(Model.numComp,1)); + % mrdm.comp = repmat(Model.constantParams{2}',2,1); + % mrdm.name = {A.name{:};B.name{:}}; + % mrdm.region = repmat(r,size(mrdm.set)); + % mrdm.hemis = repmat(h,size(mrdm.set)); + % mRDM = addstruct(mRDM,mrdm); + % % + + D.maxidx = repmat(m,length(indx),1); + D.omega = omega; + D.omegaN = omegaN; + D.logtheta = logtheta; + D.logEvidence = logEvidence; + %D.logEvidenceSplit = logEvidenceSplit; + + TT = addstruct(TT,D); + + %sh1_rsa('plot_omega_prior',D); + %sh1_rsa('plot_prior',D); + %sh1_rsa('plot_tau',D); + end + figure; + subplot(3,2,[1 2]) + lineplot(TT.maxidx,TT.logEvidence,'style_shade');ylabel('logEvidence') + + subplot(3,2,[3 4]); + barplot(TT.maxidx,TT.logtheta(:,5:end),... + 'subset',ismember(TT.maxidx,[1:10:100]));ylabel('logTau') + + subplot(3,2,[5 6]); + barplot(TT.maxidx,TT.omegaN,... + 'subset',ismember(TT.maxidx,[1:10:100]));ylabel('Omega') + + + varargout={TT,mRDM}; case 'fit_model_EB_nonlin_singlemodel' % fit single model to compare log-evidence glm = 1; useolddata = 0; @@ -537,9 +558,10 @@ sh1_rsa('plot_tau',T); varargout={T}; case 'fit_model_EB_nonlin_multiguess' - glm = 1; - useolddata = 1; - modelTerms = [1 2 4 6]; + glm = 1; + useolddata = 1; + modelTerms = [1 2 4 6]; + Niter = 100; vararginoptions(varargin(:),{'glm','useolddata','modelTerms'}) % load data set @@ -553,7 +575,7 @@ end %regions=unique(T.region); - regions=[1:8]'; + regions=[1:2]'; mRDM = []; for r=regions' for h=[1 2] @@ -563,7 +585,7 @@ Sigma(:,:,s)=reshape(D.Sigma(s,:),8,8); end; - for iter = 1:50 + for iter = 1:Niter iter initialguess = rand(size(modelTerms))*100-50; initialguess(end) = 0; @@ -593,7 +615,9 @@ % normalize omega with mean of each RDM M = feval(Model.fcn,lt(Model.numPrior+1:Model.numPrior+Model.numNonlin),... Model.constantParams{:}); - omega(indx,:) = w.*permute(mean(M.RDM,2),[3 1 2]); + normX = sqrt(mean(M.RDM.^2,2)); + omega(indx,:) = w; + omegaN(indx,:) = bsxfun(@times,w,normX'); % save resultant RDMs with estimated tau A = feval(Model.fcn,lt(Model.numPrior+1:Model.numPrior+Model.numNonlin),... @@ -694,6 +718,47 @@ imagesc_rectangle([mean(Ldata.logtheta(:,1:4),1),mean(Rdata.logtheta(:,1:4),1)],... 'scale',[-scale, scale]); end + case 'plot_reliability_btw' % plot between subject reliabitliy of RDM + C = varargin{1}; + region = [1:8]; + barwidth = 0.8; + + Nregion = numel(region); + figure('name',sprintf('Between subject reliability of RDM'),'color','w'); + for reg = region; + count=0; + for h=[2 1] + % right hemi + subplot(2,Nregion,reg+count*Nregion); + subset = C.region==reg&C.hemis==h; + upperbound = nanmean(C.within_u(subset,:)); + se = nanmean(C.within_u(subset,:))/sqrt(nancount(C.within_u(subset,:))); + + xpos = barplot([],[C.within_l,C.across,C.within_l_reduced],'subset',subset,'barwidth',barwidth); + [tval, p]=ttest_mc(C.within_l,C.across,1,'paired','subset',subset) + + drawline(0,'dir','horz','lim',[xpos(1)-barwidth,xpos(2)+barwidth]); hold on + drawline(upperbound,'dir','horz','color',[0.5 0.5 0.5],'lim',... + [xpos(1)-barwidth,xpos(2)+barwidth],... + 'linewidth',2,'error',se,'errorcolor',[0.8 0.8 0.8]) + xpos = barplot([],[C.within_l,C.across,C.within_l_reduced],'subset',subset,'barwidth',barwidth); hold on + + title([regname{reg},'-',hemName{h}],'fontsize',12); + ylabel('Pearson r','fontsize',12);xlabel(''); + set(gca,'ylim',[-0.2 0.8],'xticklabel',{'W','B'},'fontsize',12) + if p<0.05 + drawline(0.78,'dir','horz','lim',xpos([1,2]),'color',[1 0 0]); + elseif p<0.01 + drawline(0.78,'dir','horz','lim',xpos([1,2]),'color',[0 1 0]); + elseif p<0.001 + drawline(0.78,'dir','horz','lim',xpos([1,2]),'color',[0 0 1]); + end + count=count+1; + end + end + + + case 'showRDMs' T=varargin{1}; From c76319342b8121e3b288929fdb09cf0481cd56f3 Mon Sep 17 00:00:00 2001 From: Atsushi Yokoi Date: Wed, 7 Oct 2015 17:24:12 -0400 Subject: [PATCH 3/3] Changes on normalisation of omega. --- sh1_rsa.m | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/sh1_rsa.m b/sh1_rsa.m index a8445b4..4ba2096 100644 --- a/sh1_rsa.m +++ b/sh1_rsa.m @@ -386,9 +386,9 @@ glm = 3; useolddata = 0; modelTerms = [1 2 4 6]; - Niter = 100; % iteration from different initial parameter values - Ntake = 100; % how many best models you choose - region = [1 2]; + Niter = 50; % iteration from different initial parameter values + Ntake = 50; % how many best models you choose + region = [3 4]; vararginoptions(varargin(:),{'glm','useolddata','modelTerms','Niter','Ntake','region'}) TT=[]; @@ -447,7 +447,7 @@ Model.constantParams{:}); normX = sqrt(mean(M.RDM.^2,2)); omega(:,:) = ww; - omegaN(:,:) = bsxfun(@times,ww,normX'); + omegaN(:,:) = ww.*permute(normX,[3 1 2]);%bsxfun(@times,ww,normX'); % save resultant RDMs with estimated tau % A = feval(Model.fcn,lt(Model.numPrior+1:Model.numPrior+Model.numNonlin),... @@ -491,10 +491,10 @@ varargout={TT,mRDM}; case 'fit_model_EB_nonlin_singlemodel' % fit single model to compare log-evidence - glm = 1; + glm = 3; useolddata = 0; modelTerms = [1 2 4 6]; - nonlinP0 = [-10 1 10 0]; + nonlinP0 = [0 0 0 0]; vararginoptions(varargin(:),{'glm','useolddata','modelTerms','nonlinP0'}) % load data set @@ -508,7 +508,7 @@ end %regions=unique(T.region); - regions=[1:8]'; + regions=[2 3 4 5 7]'; for r=regions' fprintf('Estimating ROI distance %s...\n',regname{r}); for h=[1 2] @@ -535,12 +535,12 @@ logprior(indx,m) = repmat(lt(1),length(indx),1); logtau(indx,m) = repmat(lt(2),length(indx),1); - term(indx,m) = repmat(modelTerms(m),length(indx),1);; + term(indx,m) = repmat(modelTerms(m),length(indx),1); % normalize omega with mean of each RDM M = feval(Model.fcn,lt(Model.numPrior+1:Model.numPrior+Model.numNonlin),... Model.constantParams{:}); - omega(indx,m) = w.*permute(mean(M.RDM,2),[3 1 2]); + omega(indx,m) = w.*permute(sqrt(mean(M.RDM.^2,2)),[3 1 2]); end % model term end; @@ -726,6 +726,7 @@ Nregion = numel(region); figure('name',sprintf('Between subject reliability of RDM'),'color','w'); for reg = region; + fprintf('%s',regname{reg}) count=0; for h=[2 1] % right hemi