Unverified Commit f2866530 authored by mayucel's avatar mayucel Committed by GitHub
Browse files

GLM tstats, hmrS_GLM (#59) 

* Update SubjClass.m

dc, Aaux, tIncAuto and rcmap are also collected for each run for the subject level GLM

* Update hmrR_GLM.m

Adding tstats per condition and per contrast for the iterative restricted maximum likelihood method from Barker, BOE, 2013.

* Create hmrS_GLM.m

Adding hmrS_GLM which gets all run level data from a single subject and concatenates them. The new design matrix has HRF regressors common to all runs but have separate regressors for the rest of the model (e.g. drift, short separation).
parent 10832b7e

DataTree/SubjClass.m

+6 −0
Original line number Diff line number Diff line
@@ -249,6 +249,12 @@ classdef SubjClass < TreeNodeClass 
            obj.outputVars.mlActRuns{r.iRun}     = r.procStream.output.GetVar('mlActAuto');

            obj.outputVars.nTrialsRuns{r.iRun}   = r.procStream.output.GetVar('nTrials');

            obj.outputVars.stimRuns{r.iRun}      = r.GetVar('stim');

            obj.outputVars.dcRuns{r.iRun}       = r.procStream.output.GetVar('dc');

            obj.outputVars.AauxRuns{r.iRun}      = r.procStream.output.GetVar('Aaux');

            obj.outputVars.tIncAutoRuns{r.iRun}  = r.procStream.output.GetVar('tIncAuto');

            obj.outputVars.rcMapRuns{r.iRun}     = r.procStream.output.GetVar('rcMap');



            

            

            % a) Find all variables needed by proc stream

            args = obj.procStream.GetInputArgs();

FuncRegistry/UserFunctions/hmrR_GLM.m

+36 −20
Original line number Diff line number Diff line
@@ -608,17 +608,15 @@ for iBlk=1:length(data_y) 
                    ytmp = y(lstInc,conc,lstML);

                    for chanIdx=1:length(lstML)

                        ytmp2 = y(lstInc,conc,lstML(chanIdx));

                        [dmoco, beta, tstat, pval0, sigma, CovB, dfe, w, P, f] = ar_glm_final(squeeze(ytmp2),At(lstInc,:));

                        [dmoco, beta, tstat(:,lstML(chanIdx),conc), pval(:,lstML(chanIdx),conc), sigma, CovB(:,:,lstML(chanIdx),conc), dfe, w, P, f] = ar_glm_final(squeeze(ytmp2),At(lstInc,:));



                        foo(:,lstML(chanIdx),conc)=beta;

                        ytmp(:,1,chanIdx) = dmoco; %We also need to keep my version of "Yvar" and "Bvar"                    

                        

                        yvar(:,lstML(chanIdx),conc)=sigma.^2;

                        bvar(:,lstML(chanIdx),conc)=diag(CovB);  %Note-  I am only keeping the diag terms.  This lets you test if beta != 0,

                        bvar(:,lstML(chanIdx),conc)=diag(CovB(:,:,lstML(chanIdx),conc));  %Note-  I am only keeping the diag terms.  This lets you test if beta != 0,

                        %but in the future the HOMER-2 code needs to be modified to keep the entire cov-beta matrix which you need to test between conditions e.g. if beta(1) ~= beta(2)

                   

                        % GLM stats for each condition

                        tval(:,lstML(chanIdx),conc) = tstat;

                        pval(:,lstML(chanIdx),conc) = pval0;

                        % reply to above comment from DAB: Now directly using CovB for the contrast calculations. MAY

                    end

                end
@@ -685,10 +683,28 @@ for iBlk=1:length(data_y) 
                        end

                    end

                    

                elseif glmSolveMethod==2  % WLS

                elseif glmSolveMethod==2  % iWLS

                    

                    yest(:,lstML,conc) = At * foo(:,lstML,conc);

                    for iCh = 1:length(lstML)

                        

                        % GLM stats for contrast between conditions, given a c_vector exists

                        if nCond > 1

                            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(At,2)-size(cv_dummy,2))];

                                end

                                tval_contrast(:,lstML(iCh),conc) = cv_extended * foo(:,lstML(iCh),conc)./sqrt(cv_extended * squeeze(CovB(:,:,lstML(chanIdx),conc))* cv_extended');

                                pval_contrast(:,lstML(iCh),conc) = 1-tcdf(abs(tval_contrast(:,lstML(iCh),conc)),(size(y,1)-1));

                            end

                        end

                        %

                        

                        for iCond=1:nCond

                            if size(tbasis,3)==1

                                yavgstd(:,lstML(iCh),conc,iCond) = diag(tbasis*diag(bvar([1:nB]+(iCond-1)*nB,lstML(iCh),conc))*tbasis').^0.5;
@@ -785,24 +801,24 @@ for iBlk=1:length(data_y) 
    yR_blks{iBlk}   = yR;

    

 % stats struct

    if glmSolveMethod == 1 %  for OLS only for now

    if glmSolveMethod == 1 % for OLS

        % GLM stats for each condition

        hmrstats.beta_label = beta_label;

        hmrstats.tval = tval;

        hmrstats.pval = pval;

        hmrstats.ml = ml;

    else                   % for iWLS

        hmrstats.beta_label = beta_label;

        hmrstats.tval = tstat;

        hmrstats.pval = pval;

        hmrstats.ml = ml;

    end

    

    % GLM stats for contrast between conditions, if c_vector exists

    if (sum(abs(c_vector)) ~= 0) && (size(c_vector,2) == nCond) && nCond>1

        hmrstats.tval_contrast = tval_contrast;

        hmrstats.pval_contrast = pval_contrast;

        hmrstats.contrast = c_vector;

    end

    elseif glmSolveMethod == 2 %  for AR 

                % GLM stats for each condition 

        hmrstats.beta_label = beta_label;

        hmrstats.tval = tval;

        hmrstats.pval = pval;

        hmrstats.ml = ml;        

    end

    

end

FuncRegistry/UserFunctions/hmrS_GLM.m

0 → 100644
+900 −0

File added.

Preview size limit exceeded, changes collapsed.

CRA Git | Maintained and supported by SUSTech CRA and CCSE