Update hmrR_GLM.m (014edcaa) · Commits · github_fork / Homer3

FuncRegistry/UserFunctions/hmrR_GLM.m

+77 −2

Original line number	Diff line number	Diff line
		@@ -354,6 +354,7 @@ for iBlk=1:length(data_y)
		end
		clmn = clmn(1:nT);
		dA(:,iC,iConc) = clmn;
		beta_label{b + (iCond-1)*nB} = ['Cond' num2str(iCond)];
		end
		end
		end
		@@ -403,14 +404,42 @@ for iBlk=1:length(data_y)
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		% Final design matrix
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		dummy = size(beta_label,2);

		switch flagNuisanceRMethod
		case {0,1,2} % short separation
		for iConc=1:2
		A(:,:,iConc)=[dA(:,:,iConc) xDrift Aaux Amotion];

		if iConc == 1
		for ixDrift = 1:size(xDrift,2)
		beta_label{ixDrift + dummy} = ['xDrift'];
		end
		dummy = size(beta_label,2);
		for iAaux = 1:size(Aaux,2)
		beta_label{iAaux + dummy} = ['Aux'];
		end
		dummy = size(beta_label,2);
		for iAmotion = 1:size(Amotion,2)
		beta_label{ixDrift + dummy} = ['Motion'];
		end
		end

		end
		case 3 % tCCA regressor design matrix without Aaux (will be put in in the loop below)
		for iConc=1:2
		A(:,:,iConc)=[dA(:,:,iConc) xDrift Amotion];

		if iConc == 1
		for ixDrift = 1:size(xDrift,2)
		beta_label{ixDrift + dummy} = ['Drift'];
		end
		dummy = size(beta_label,2);
		for iAmotion = 1:size(Amotion,2)
		beta_label{ixDrift + dummy} = ['Motion'];
		end
		end

		end
		end

		@@ -468,24 +497,64 @@ for iBlk=1:length(data_y)
		% lstML = find(iNearestSS==mlSSlst(iSS));
		Ass = y(:,conc,mlSSlst(iSS));
		At = [A(:,:,conc) Ass];

		if iSS == 1 && conc == 1
		dummy = size(beta_label,2);
		for iAss = 1:size(Ass,2)
		beta_label{iAss + dummy} = ['ShortSep'];
		end
		end

		elseif flagNuisanceRMethod==1
		lstML = find(iNearestSS(:,conc)==mlSSlst(iSS) & mlAct(lstMLtmp)==1);
		% lstML = find(iNearestSS(:,conc)==mlSSlst(iSS));
		Ass = y(:,conc,mlSSlst(iSS));
		At = [A(:,:,conc) Ass];

		if iSS == 1 && conc == 1
		dummy = size(beta_label,2);
		for iAss = 1:size(Ass,2)
		beta_label{iAss + dummy} = ['ShortSep'];
		end
		end

		elseif flagNuisanceRMethod==2
		lstML = lstMLtmp(find(mlAct(lstMLtmp)==1));
		% lstML = 1:size(y,3);
		Ass = mean(y(:,conc,lstSS),3);
		At = [A(:,:,conc) Ass];

		if iSS == 1 && conc == 1
		dummy = size(beta_label,2);
		for iAss = 1:size(Ass,2)
		beta_label{iAss + dummy} = ['ShortSep'];
		end
		end

		elseif flagNuisanceRMethod==3
		if ischar(rcMap{iBlk}) % no channel map: use all tCCA regressors for one group of all channels
		lstML = lstMLtmp(find(mlAct(lstMLtmp)==1));
		At = [A(:,:,conc) Aaux];

		if iSS == 1 && conc == 1
		dummy = size(beta_label,2);
		for iAaux = 1:size(Aaux,2)
		beta_label{iAaux + dummy} = ['Aux'];
		end
		end

		elseif iscell(rcMap{iBlk}) % channel map: each single regressor corresponds to one channel (nCH groups)
		lstML = lstMLtmp(find(mlAct(rcMap{iBlk}{conc,iSS})==1));
		Atcca = Aaux{iBlk}(:,rcMap{conc,iSS});
		At = [A(:,:,conc) Atcca];

		if iSS == 1 && conc == 1
		dummy = size(beta_label,2);
		for iAtcca = 1:size(Atcca,2)
		beta_label{iAtcca + dummy} = ['tCCA regressor'];
		end
		end

		end
		end

		@@ -573,11 +642,12 @@ for iBlk=1:length(data_y)
		if glmSolveMethod==1 % OLS ~flagUseTed
		pAinvAinvD = diag(pinvA*pinvA');
		yest(:,lstML,conc) = At * foo(:,lstML,conc);
		yvar(1,lstML,conc) = std(squeeze(y(:,conc,lstML))-yest(:,lstML,conc),[],1).^2; % check this against eq(53) in Ye2009
		yvar(1,lstML,conc) = sum((squeeze(y(:,conc,lstML))-yest(:,lstML,conc)).^2)./(size(y,1)-1); % check this against eq(53) in Ye2009
		for iCh = 1:length(lstML)
		bvar(:,lstML(iCh),conc) = yvar(1,lstML(iCh),conc) * pAinvAinvD;
		tval(:,lstML(iCh),conc) = foo(:,lstML(iCh),conc)./sqrt(bvar(:,lstML(iCh),conc));
		pval(:,lstML(iCh),conc) = 1-tcdf(abs(tval(:,lstML(iCh),conc)),(size(y,1)-1));
		for iCond=1:nCond
		% yavgstd(:,lstML(iCh),conc,iCond) = diag(tbasisdiag(bvar([1:nB]+(iCond-1)nB,lstML(iCh),conc))*tbasis').^0.5;
		if size(tbasis,3)==1
		yavgstd(:,lstML(iCh),conc,iCond) = diag(tbasisdiag(bvar([1:nB]+(iCond-1)nB,lstML(iCh),conc))*tbasis').^0.5;
		else
		@@ -686,5 +756,10 @@ for iBlk=1:length(data_y)
		beta_blks{iBlk} = beta;
		yR_blks{iBlk} = yR;

		% stats struct
		hmrstats.beta_label = beta_label; %
		hmrstats.tval = tval;
		hmrstats.pval = pval;

		end

Admin message