54. Genearal linear model on fMRI data (individual and group analysis) - fahsuanlin/labmanual GitHub Wiki

These are procedures modeling the fMRI time series using General Linear Model (GLM) with more detailed scripts in analysis 1-st level individual results and 2nd-level group results. Consult the introductory fMRI analysis page. Data were tested by the Human Connectome Project data. Since this is task-fMRI analysis, we will take the working memory (WM) data as the example.

NOTE: At Sunnybrook, the data and scripts are at /space_lin1/fhlin/hcp and /space_lin1/fhlin/hcp/scripts, respectively.

1. 1-st level analysis on an individual

Following the described procedures for GLM analysis.

1.1 Calculate the effect and statistics related to a single hypothesis.

Consider the example script here. Use t-statistics and save the associated effect as the following:

• Calculate the coefficients for regressors using least squares.

        fprintf('GLM estimation....');
        beta=inv(contrast'*contrast)*contrast'*stc;
        error=(stc-contrast*beta);
        error_sig2=sum(error.^2,1)./(size(stc,1)-size(contrast,2));
        df=size(stc,1)-rank(contrast);

• Calculate the effect and the associated t-statistics

            cc=zeros(size(w_beta,1),1);
            cc_tmp=kron(cc_hdr(:),ones(size(HDR0,2),1));
            cc(1:length(cc_tmp(:)))=cc_tmp(:);
            t_stat=cc'*w_beta./sqrt(w_error_sig2.*(cc'*inv(w_contrast'*w_contrast)*cc));
            effect=cc'*w_beta;

NOTE: w_beta, w_error_sig2, and w_contrast are identical to beta, error_sig2, and contrast, respectively, if no auto-regressive modeling on the fMRI time series.

1.2 Calculate the effect and statistics related to multiple hypothesis.

Consider the example script here. Use F-statistics and save the associated effect as the following:

• Calculate the coefficients for regressors using least squares.

        fprintf('GLM estimation....');
        beta=inv(contrast'*contrast)*contrast'*stc;
        error=(stc-contrast*beta);
        error_sig2=sum(error.^2,1)./(size(stc,1)-size(contrast,2));
        df=size(stc,1)-rank(contrast);

• Calculate the effect and the associated F-statistics. Here effects are multiple contrast images, each of which corresponds to one HRF basis.

            %effect from F-test (more than 1 HRF basis)
            cc_tmp=kron(cc_hdr(:),eye(size(HDR0,2)))';
            for cc_idx=1:size(cc_tmp,1)
                cc=zeros(size(w_beta,1),1);
                cc(1:size(cc_tmp,2))=cc_tmp(cc_idx,:);
                t_stat(cc_idx,:)=cc'*w_beta./sqrt(w_error_sig2.*(cc'*inv(w_contrast'*w_contrast)*cc));
                effect(cc_idx,:)=cc'*w_beta;
            end;

NOTE: w_beta, w_error_sig2, and w_contrast are identical to beta, error_sig2, and contrast, respectively, if no auto-regressive modeling on the fMRI time series.

2. 2-nd level analysis on an individual

2.1 Calculate the statistics related for a single hypothesis.

Consider the example script here. The following are codes to get Z-scores for a single hypothesis. Basically it's just averaging.

               for d_idx=1:length(subject)
                        fprintf('{%s}::[%s]...(%04d|%04d)....\r',file_stem{f_idx},subject{d_idx},d_idx,length(subject));

                        fn=sprintf('%s/%s/analysis/%s-%s.stc',data_path,subject{d_idx},file_stem{f_idx},hemi_str);
                        if(exist(fn))
                                [stc,v]=inverse_read_stc(fn);

                                stc_all(f_idx).stc(:,counter,hemi)=stc(:,1);
                                if(counter==1)
                                        tmp=stc(:,1);
                                        tmp2=stc(:,1).^2;
                                else
                                        tmp=tmp+stc(:,1);
                                        tmp2=tmp2+stc(:,1).^2;
                                end;
                                counter=counter+1;
                        end;
                end;
                fprintf('\n');
                stc_avg=tmp./counter;
                tmp2=tmp2./counter;
                stc_std=sqrt(tmp2-stc_avg.^2);
                stc_z=stc_avg./stc_std.*sqrt(counter-1);

You can also implement another GLM on saved individual effects as the 1-st level analysis. See 1.1 above.

2.2 Calculate the statistics related to a composite hypothesis.

Consider the example script here. An F-statistics is calculated.

        for v_idx=1:size(stc_all,1)
            buffer=squeeze(stc_all(v_idx,:,:))';


            D=ones(size(buffer,1),1); %average....

            beta=inv(D'*D)*D'*buffer;
            res=buffer-D*beta;


	    S=res'*res./size(res,1);

	    T2=beta*inv(S)*beta'.*size(res,1); %hotelling's t squares
	    F(v_idx)=T2.*(size(res,1)-size(beta,2))./(size(res,1)-1)./size(beta,2);
        end;
        dof_num=size(buffer,2);
        dof_den=size(res,1)-size(buffer,2);

stc_all is a 3D data matrix of effects across individuals (#voxel x #hypothesis x #individual). In this example, we test a composite null hypothesis of all hypothesis effect is zero. The resultant F is the group-level statistics with degrees of freedom dof_num and dof_den.