% Plotting the hemodynamic response function (hrf) using fmridesign hrf_parameters=[5.4 5.2 10.8 7.35 0.35] time=(0:240)/10; hrf0=fmridesign(time,0,[1 0],[],hrf_parameters); plot(time,squeeze(hrf0(:,1,1)),'LineWidth',2) xlabel('time (seconds)') ylabel('hrf') title('Glover hrf model') % Making the design matrices using fmridesign frametimes=(0:119)*3; slicetimes=[0.14 0.98 0.26 1.10 0.38 1.22 0.50 1.34 0.62 1.46 0.74 1.58 0.86]; eventid=kron(ones(10,1),[1; 2]); eventimes=(0:19)'*18+9; duration=ones(20,1)*9; height=ones(20,1); events=[eventid eventimes duration height] S=kron(ones(10,1),kron([0 0; 1 0; 0 0; 0 1],ones(3,1))); X_cache=fmridesign(frametimes,slicetimes,events,[],hrf_parameters); X_cache=fmridesign(frametimes,slicetimes, [] , S,hrf_parameters); plot(squeeze(X_cache(:,:,1,5)),'LineWidth',2) legend('Hot','Warm') xlabel('frame number') ylabel('response') title('Hot and Warm responses') % Analysing one run with fmrilm contrast=[1 0; 0 1; 1 -1]; exclude=[1 2]; which_stats=[1 1 1 1 1]; input_file='e:/keith/data/brian_ha_19971124_1_100326_mri_MC.mnc'; output_file_base=['e:/keith/results/ha_100326_hot'; 'e:/keith/results/ha_100326_wrm'; 'e:/keith/results/ha_100326_hmw']; [df1 p]=fmrilm(input_file, output_file_base, X_cache, contrast, exclude, which_stats); m=fmris_read_image('e:/keith/results/ha_100326_hmw_mag_t.mnc',4,1); imagesc(m.data'); colorbar; axis xy; a=fmris_read_minc(input_file,4,1); imagesc((m.data.*(a.data>500))'); colorbar; axis xy; m=fmris_read_image('e:/keith/results/ha_100326_hmw_mag_t.mnc',4:6,1); mm=reshape(permute(m.data,[2 1 3 4]),m.dim(1),prod(size(m.data))/m.dim(1)); imagesc(mm.*repmat(a.data'>500,1,3)); colorbar; axis xy; contrast=[0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 1] which_stats=[0 0 0 1]; output_file_base='e:/keith/results/ha_100326_hot'; fwhm_rho='e:/keith/results/ha_100326_hot_rho.mnc'; fmrilm(input_file,output_file_base,X_cache,contrast,exclude,which_stats,fwhm_rho); % A linear effect of temperature temperature=[45.5 35.0 49.0 38.5 42.0 49.0 35.0 42.0 38.5 45.5 ... 38.5 49.0 35.0 45.5 42.0 45.5 38.5 42.0 35.0 49.0]'; Temperature=35:3.5:49; subplot(2,2,1) plot(Temperature,1:5,'LineWidth',2); xlabel('Temperature'); ylabel('Response'); title('(a) Linear temperature effect'); subplot(2,2,2) plot(Temperature,10-((1:5)-4).^2,'LineWidth',2); xlabel('Temperature'); ylabel('Response'); title('(b) Quadratic temperature effect'); subplot(2,2,3) plot(Temperature,10-((1:5)-4).^2+((1:5)-3).^3,'LineWidth',2); xlabel('Temperature'); ylabel('Response'); title('(c) Cubic temperature effect'); subplot(2,2,4) plot(Temperature,[1 4 3 5 3.5],'LineWidth',2); xlabel('Temperature'); ylabel('Response'); title(['(d) Quartic = arbitrary temperature effect']); events=[zeros(20,1)+1 eventimes duration ones(20,1); zeros(20,1)+2 eventimes duration temperature] contrast=[0 1; 1 49; 1 35; 0 14]; events=[zeros(20,1)+1 eventimes duration ones(20,1); zeros(20,1)+2 eventimes duration temperature; zeros(20,1)+3 eventimes duration temperature^2] contrast=[0 0 1]; contrast=[1 0 0 0 0; 0 1 0 0 0; 0 0 1 0 0; 0 0 0 1 0; 0 0 0 0 1]; contrast=[eye(5) ones(5,4)]; which_stats=[0 0 0 1]; contrast=[.8 -.2 -.2 -.2 -.2; -.2 .8 -.2 -.2 -.2; -.2 -.2 .8 -.2 -.2; -.2 -.2 -.2 .8 -.2; -.2 -.2 -.2 -.2 .8]; contrast=[eye(5)-ones(5)/5 ones(5,4)]; which_stats=[0 0 0 1]; contrast=[-7.0 -3.5 0 3.5 7.0]; which_stats=[1 1 1]; % Combining runs/sessions/subjects with multistat contrast=[1 -1]; which_stats=[1 1 1]; input_file='e:/keith/data/brian_ha_19971124_1_093923_mri_MC.mnc'; output_file_base='e:/keith/results/ha_093923_hmw'; [df2 p]=fmrilm(input_file, output_file_base, X_cache, contrast, exclude, which_stats); input_file='e:/keith/data/brian_ha_19971124_1_101410_mri_MC.mnc'; output_file_base='e:/keith/results/ha_101410_hmw'; [df3 p]=fmrilm(input_file, output_file_base, X_cache, contrast, exclude, which_stats); input_file='e:/keith/data/brian_ha_19971124_1_102703_mri_MC.mnc'; output_file_base='e:/keith/results/ha_102703_hmw'; [df4 p]=fmrilm(input_file, output_file_base, X_cache, contrast, exclude, which_stats); X=[1 1 1 1]'; contrast=[1]; which_stats=[1 1 1]; input_files_effect=['e:/keith/results/ha_093923_hmw_mag_ef.mnc'; 'e:/keith/results/ha_100326_hmw_mag_ef.mnc'; 'e:/keith/results/ha_101410_hmw_mag_ef.mnc'; 'e:/keith/results/ha_102703_hmw_mag_ef.mnc']; input_files_sdeffect=['e:/keith/results/ha_093923_hmw_mag_sd.mnc'; 'e:/keith/results/ha_100326_hmw_mag_sd.mnc'; 'e:/keith/results/ha_101410_hmw_mag_sd.mnc'; 'e:/keith/results/ha_102703_hmw_mag_sd.mnc']; output_file_base='e:/keith/results/ha_multi_hmw' input_files_df=[df1 df2 df3 df4] input_files_fwhm=6 df=multistat(input_files_effect,input_files_sdeffect,input_files_df, ... input_files_fwhm,X,contrast,output_file_base,which_stats,Inf) % Fixed and random effects which_stats=[1 1 1 1]; df=multistat(input_files_effect,input_files_sdeffect,input_files_df, ... input_files_fwhm,X,contrast,output_file_base,which_stats,15) m=fmris_read_image('e:/keith/results/ha_multi_hmw_sdratio.mnc',4,1); imagesc(m.data'); colorbar; axis xy; % Thresholding the tstat image with stat_threshold and fdr_threshold stat_threshold(1000000,38.4521,6,112,0.05) stat_threshold(1000000,38.4521,6,112,0.05,3) mask_file='e:/keith/data/brian_ha_19971124_1_100326_mri_MC.mnc'; input_file='e:/keith/results1/ha_multi_hmw_t.mnc'; fdr_threshold(input_file,[], mask_file,500,112,0.05) % Locating peaks with locmax lm=locmax('e:/keith/results/ha_multi_hmw_t.mnc',4.86) pval=stat_threshold(1000000,38.4521,6,112,lm(:,1)) % Extracting values from a minc file using extract voxel=[62 65 4] effect=extract(voxel,'e:/keith/results/ha_multi_hmw_ef.mnc') sdeffect=extract(voxel,'e:/keith/results/ha_multi_hmw_sd.mnc') ref_data=squeeze(extract(voxel,'e:/keith/data/brian_ha_19971124_1_100326_mri_MC.mnc')); ef_hot=extract(voxel,'e:/keith/results/ha_100326_hot_mag_ef.mnc') ef_wrm=extract(voxel,'e:/keith/results/ha_100326_wrm_mag_ef.mnc') fitted=mean(ref_data)+ef_hot*X_cache(:,1,1,5)+ef_wrm*X_cache(:,2,1,5); plot(frametimes,[ref_data fitted],'LineWidth',2); ylim([920 960]); legend('Reference data','Fitted values'); xlabel('time (seconds) from start of epoch'); ylabel('fMRI response'); % Estimating the time course of the response eventid=kron(ones(10,1),(1:12)'); eventimes=frametimes'; duration=ones(120,1)*3; height=ones(120,1); events=[eventid eventimes duration height] X_hrf=fmridesign(frametimes,slicetimes,events,[],zeros(1,5)); plot(squeeze(X_hrf(:,1:4,1,5)),'LineWidth',2) legend('Time01','Time02','Time03','Time04') xlabel('frame number') ylabel('response') title('Slice 5, first 4 times only') contrast=[eye(12)-ones(12)/12]; exclude=[1 2]; which_stats=[0 1 1 1]; input_file='e:/keith/data/brian_ha_19971124_1_100326_mri_MC.mnc'; output_file_base=['e:/keith/results/ha_100326_time01'; 'e:/keith/results/ha_100326_time02'; 'e:/keith/results/ha_100326_time03'; 'e:/keith/results/ha_100326_time04'; 'e:/keith/results/ha_100326_time05'; 'e:/keith/results/ha_100326_time06'; 'e:/keith/results/ha_100326_time07'; 'e:/keith/results/ha_100326_time08'; 'e:/keith/results/ha_100326_time09'; 'e:/keith/results/ha_100326_time10'; 'e:/keith/results/ha_100326_time11'; 'e:/keith/results/ha_100326_time12']; fmrilm(input_file,output_file_base,X_cache,contrast,exclude,which_stats,fwhm_rho); stat_threshold(1000000,38.4538,6,[11 103]) lm=locmax([output_file_base(1,:) '_mag_F.mnc'],5.18) voxel=[63 66 4] values=extract(voxel,output_file_base,'_mag_ef.mnc') sd=extract(voxel,output_file_base,'_mag_sd.mnc') b_hot=extract(voxel,'e:/keith/results/ha_100326_hot_mag_ef.mnc') b_wrm=extract(voxel,'e:/keith/results/ha_100326_wrm_mag_ef.mnc') time=(1:360)/10; X_hrf=fmridesign(time,0,[1 9 9 1],[],hrf_parameters); hrf=squeeze(X_hrf(:,1,1)); plot((0:12)*3+slicetimes(voxel(3)),values([1:12 1]),'k', ... [0:11; 0:11]*3+slicetimes(voxel(3)), [values+sd; values-sd],'g', ... time,[zeros(1,90) ones(1,90) zeros(1,180)],'r', ... time,[zeros(1,270) ones(1,90)],'b', ... time,hrf*b_hot+hrf([181:360 1:180])*b_wrm,'g','LineWidth',2); legend('Estimated response','Modeled response'); text(10,0.5,'Hot') text(28,0.5,'Warm') xlabel('time (seconds) from start of epoch'); ylabel('fMRI response'); % Estimating the delay contrast=[1 0; 0 1; 1 -1] which_stats=[1 1 1] input_file='e:/keith/data/brian_ha_19971124_1_100326_mri_MC.mnc'; output_file_base=['e:/keith/results/ha_100326_hot'; 'e:/keith/results/ha_100326_wrm'; 'e:/keith/results/ha_100326_hmw'] n_poly=3; fwhm_rho='e:/keith/results/ha_100326_hot_rho.mnc'; reference_delay=5.4; contrast_is_delay=[1 1 1]; fmrilm(input_file, output_file_base, X_cache, contrast, exclude, ... which_stats, fwhm_rho, n_poly, reference_delay, contrast_is_delay); t=fmris_read_minc('e:/keith/results/ha_100326_hot_mag_t.mnc',4,1); d=fmris_read_minc('e:/keith/results/ha_100326_hot_del_ef.mnc',4,1); s=fmris_read_minc('e:/keith/results/ha_100326_hot_del_sd.mnc',4,1); a=fmris_read_minc(input_file,4,1); imagesc([(d.data.*(t.data>4)+(a.data<500))' (s.data.*(t.data>4)+(a.data<500))']); colorbar; axis xy; % Effective connectivity of all voxels with a reference voxel input_file='e:/keith/data/brian_ha_19971124_1_100326_mri_MC.mnc'; output_file_base='e:/keith/results/ha_100326_connect'; contrast=[0 0 0 0 0 0 1]; which_stats=[1 1 1]; reference_delay=0; contrast_is_delay=0; confounds=ref_data; fmrilm(input_file, output_file_base, X_cache, contrast, exclude, ... which_stats, fwhm_rho, n_poly, reference_delay, contrast_is_delay, confounds, 2); m=fmris_read_image('e:/keith/results/ha_100326_connect_mag_t.mnc',4,1); imagesc((m.data.*(a.data>500))'); colorbar; axis xy; % Higher order autoregressive models which_stats=[1 1 1 0 1 0 0 1]; fmrilm(input_file, output_file_base, X_cache, contrast, exclude, ... which_stats, 15, 3, 0, 0, [], 2); m=fmris_read_minc('e:/keith/results/ha_100326_hot_A.mnc',4,1:2); mm=reshape(permute(m.data,[2 1 3 4]),m.dim(1),prod(size(m.data))/m.dim(1)); imagesc(mm); colorbar; axis xy;