minor change (053971c7) · Commits · github_fork / Homer3

FuncRegistry/UserFunctions/hmrR_GLM.m

+72 −41

Original line number	Original line	Diff line number	Diff line
	% SYNTAX:		% SYNTAX:
	% [yavg, yavgstd, tHRF, nTrials, ynew, yresid, ysum2, beta, R] = hmrR_GLM(data, stim, probe, mlActAuto, Aaux, tIncAuto, trange, rcMap, glmSolveMethod, idxBasis, paramsBasis, rhoSD_ssThresh, flagNuisanceRMethod, driftOrder)		% [yavg, yavgstd, tHRF, nTrials, ynew, yresid, ysum2, beta, R, hmrstats] = hmrR_GLM(data, stim, probe, mlActAuto, Aaux, tIncAuto, trange, rcMap, glmSolveMethod, idxBasis, paramsBasis, rhoSD_ssThresh, flagNuisanceRMethod, driftOrder, c_vector)
	%		%
	% UI NAME:		% UI NAME:
	% GLM_HRF_Drift_SS		% GLM_HRF_Drift_SS
	@@ -69,7 +69,9 @@
	% 3. uses tCCA regressors for nuisance regression, in Aaux,		% 3. uses tCCA regressors for nuisance regression, in Aaux,
	% mapped by rcMap, provided by hmr_tCCA()		% mapped by rcMap, provided by hmr_tCCA()
	% driftOrder - Polynomial drift correction of this order		% driftOrder - Polynomial drift correction of this order
	%		% c_vector - Contrast vector, has values 1, -1 or 0. E.g. to contrast cond
			% 2 to cond 3 in an experimental paradigm with four conditions, c_vector is
			% [0 1 -1 0]
	%		%
	% OUTPUTS:		% OUTPUTS:
	% yavg - the averaged results		% yavg - the averaged results
	@@ -84,10 +86,10 @@
	% (#coefficients x HbX x #Channels x #conditions)		% (#coefficients x HbX x #Channels x #conditions)
	% R - the correlation coefficient of the GLM fit to the data		% R - the correlation coefficient of the GLM fit to the data
	% (#Channels x HbX)		% (#Channels x HbX)
	%		% hmrstats - outputs t and pvalues for GLM
	%		%
	% USAGE OPTIONS:		% USAGE OPTIONS:
	% GLM_HRF_Drift_SS_Concentration: [dcAvg, dcAvgStd, nTrials, dcNew, dcResid, dcSum2, beta, R] = hmrR_GLM(dc, stim, probe, mlActAuto, Aaux, tIncAuto, rcMap, trange, glmSolveMethod, idxBasis, paramsBasis, rhoSD_ssThresh, flagNuisanceRMethod, driftOrder)		% GLM_HRF_Drift_SS_Concentration: [dcAvg, dcAvgStd, nTrials, dcNew, dcResid, dcSum2, beta, R, hmrstats] = hmrR_GLM(dc, stim, probe, mlActAuto, Aaux, tIncAuto, rcMap, trange, glmSolveMethod, idxBasis, paramsBasis, rhoSD_ssThresh, flagNuisanceRMethod, driftOrder, c_vector)
	%		%
	%		%
	% PARAMETERS:		% PARAMETERS:
	@@ -98,10 +100,10 @@
	% rhoSD_ssThresh: 15.0		% rhoSD_ssThresh: 15.0
	% flagNuisanceRMethod: 1		% flagNuisanceRMethod: 1
	% driftOrder: 3		% driftOrder: 3
			% c_vector: 0
	%		%
	%		function [data_yavg, data_yavgstd, nTrials, data_ynew, data_yresid, data_ysum2, beta_blks, yR_blks, hmrstats] = ...
	function [data_yavg, data_yavgstd, nTrials, data_ynew, data_yresid, data_ysum2, beta_blks, yR_blks] = ...		hmrR_GLM(data_y, stim, probe, mlActAuto, Aaux, tIncAuto, rcMap, trange, glmSolveMethod, idxBasis, paramsBasis, rhoSD_ssThresh, flagNuisanceRMethod, driftOrder, c_vector)
	hmrR_GLM(data_y, stim, probe, mlActAuto, Aaux, tIncAuto, rcMap, trange, glmSolveMethod, idxBasis, paramsBasis, rhoSD_ssThresh, flagNuisanceRMethod, driftOrder)

	% Init output		% Init output
	data_yavg = DataClass().empty();		data_yavg = DataClass().empty();
	@@ -112,6 +114,7 @@ data_yresid = DataClass().empty();
	beta_blks = cell(length(data_y),1);		beta_blks = cell(length(data_y),1);
	yR_blks = cell(length(data_y),1);		yR_blks = cell(length(data_y),1);
	beta_label = [];		beta_label = [];
			hmrstats = [];

	% Check input args		% Check input args
	if isempty(tIncAuto)		if isempty(tIncAuto)
	@@ -646,9 +649,30 @@ for iBlk=1:length(data_y)
	yest(:,lstML,conc) = At * foo(:,lstML,conc);		yest(:,lstML,conc) = At * foo(:,lstML,conc);
	yvar(1,lstML,conc) = sum((squeeze(y(:,conc,lstML))-yest(:,lstML,conc)).^2)./(size(y,1)-1); % 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)		for iCh = 1:length(lstML)

			% GLM stats for each condition
	bvar(:,lstML(iCh),conc) = yvar(1,lstML(iCh),conc) * pAinvAinvD;		bvar(:,lstML(iCh),conc) = yvar(1,lstML(iCh),conc) * pAinvAinvD;
	tval(:,lstML(iCh),conc) = foo(:,lstML(iCh),conc)./sqrt(bvar(:,lstML(iCh),conc));		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));		pval(:,lstML(iCh),conc) = 1-tcdf(abs(tval(:,lstML(iCh),conc)),(size(y,1)-1));
			%

			% GLM stats for contrast between conditions, given a c_vector exists
			if (sum(abs(c_vector)) ~= 0) && (size(c_vector,2) == nCond)

			if ~exist('cv_extended') == 1
			cv_dummy = [];
			for m = 1:nCond
			cv_dummy = [cv_dummy ones(1,nB)*c_vector(m)];
			end
			cv_extended = [cv_dummy zeros(1,size(beta_label,2)-size(cv_dummy,2))];
			end

			tval_contrast(:,lstML(iCh),conc) = cv_extended * foo(:,lstML(iCh),conc)./sqrt(cv_extended * (pinvApinvA') yvar(:,lstML(iCh),conc) * cv_extended');
			pval_contrast(:,lstML(iCh),conc) = 1-tcdf(abs(tval_contrast(:,lstML(iCh),conc)),(size(y,1)-1));
			end
			%


	for iCond=1:nCond		for iCond=1:nCond
	if size(tbasis,3)==1		if size(tbasis,3)==1
	yavgstd(:,lstML(iCh),conc,iCond) = diag(tbasisdiag(bvar([1:nB]+(iCond-1)nB,lstML(iCh),conc))*tbasis').^0.5;		yavgstd(:,lstML(iCh),conc,iCond) = diag(tbasisdiag(bvar([1:nB]+(iCond-1)nB,lstML(iCh),conc))*tbasis').^0.5;
	@@ -759,11 +783,18 @@ for iBlk=1:length(data_y)
	yR_blks{iBlk} = yR;		yR_blks{iBlk} = yR;

	% stats struct		% stats struct
	if glmSolveMethod == 1 % for OLS only now		if glmSolveMethod == 1 % for OLS only for now
			% GLM stats for each condition
	hmrstats.beta_label = beta_label;		hmrstats.beta_label = beta_label;
	hmrstats.tval = tval;		hmrstats.tval = tval;
	hmrstats.pval = pval;		hmrstats.pval = pval;
			hmrstats.ml = ml;
			% GLM stats for contrast between conditions, if c_vector exists
			if (sum(abs(c_vector)) ~= 0) && (size(c_vector,2) == nCond)
			hmrstats.tval_contrast = tval_contrast;
			hmrstats.pval_contrast = pval_contrast;
			hmrstats.contrast = c_vector;
	end		end

	end		end

			end
			No newline at end of file

Admin message