Dear All:
My apologies. I know that random effects analysis can be done from the GUI, but I was wondering if there were any good code examples to do this example.
See, for example, my attached code for a GLM to see what I mean. Something simple would be wonderful.
I have a bunch of contrasts images, 2 groups, and would just like to do a t test with FDR correction.
I've done it in MATLAB but no FDR correction there.
Again thank you
Sincerely yours
Misha Koshelev
p.s. Actually attached the t test MATLAB code as well just so you can see what I mean. Thank you!
--
Misha Koshelev
MD/PhD Student
Human Neuroimaging Laboratory
One Baylor Plaza
S104
Baylor College of Medicine
Houston, TX 77030
function analyze(runid,group,subj,tries)
% Analyze investors using SPM2
%
% input:
% runid - unique run identifier
% group - group identifier
% subj - subject
% tries - # of times to retry on failure
%
% Misha Koshelev
% February 24th, 2010
% Montague Laboratory
if ~exist('tries','var')
tries=3;
end
if regexp(runid,'3.*')
analyze3(runid,group,subj,tries);
return;
elseif strcmp(runid,'4')
analyze4(runid,group,subj,tries);
return;
elseif regexp(runid,'6.*')
analyze6(runid,group,subj,tries);
return;
elseif strcmp(runid,'7')
analyze7(runid,group,subj,tries);
return;
elseif strcmp(runid,'8')
analyze8(runid,group,subj,tries);
return;
end
warning off;
group
subj
% save paths
paths=path;pwds=pwd;
addpath(genpath(pwd));
succeed=false;
try
% load data
target=['/mnt/nfs/proj/trust_agnostic',filesep,group,filesep,subj];
cd([target,filesep,'other']);
load('FileInfo');
load('motion');
load('source');
source=[source,filesep,'GLMs',filesep,subj];
cd('..');
if exist(runid,'dir')
system(['rm -rf ',runid]);
end
mkdir(runid);
% image files
cd([source,filesep,'run1']);
f=spm_get('Files',pwd,'sw*.img');
cd([target,filesep,runid]);
spm_defaults;
global defaults;
% HRF
SPM.xBF.name='hrf'; %'hrf (with time derivative)';
SPM.xBF.order=1;
SPM.xBF.length=32.1250; % length in seconds
SPM.xBF.T=16; % number of time bins per scan
SPM.xBF.T0=1; % middle slice/timebin
SPM.xBF.UNITS='secs'; % OPTIONS: 'scans'|'secs' for onsets
SPM.xBF.Volterra=1; % OPTIONS: 1|2=order of convolution
% Trials
sess=1;i=1;
events=fieldnames(FileInfo.EventStruct{1});
for j=1:length(events)
event=FileInfo.EventStruct{1}.(events{j});
if strcmp(events{j},'ScreenInvAmount') % 4359 SHOWING MESSAGE <Make Initial Investment> AND PIC <0>
% rounds 1&2 only - no regressor
SPM.Sess(sess).U(i).name={[events{j},'1-2']};
SPM.Sess(sess).U(i).ons=event.values(1:2);
SPM.Sess(sess).U(i).dur=event.dur(1:2);
SPM.Sess(sess).U(i).P(1).name='none';
i=i+1;
% rounds 3-10 - regressed with i&r's at past 2 time points
SPM.Sess(sess).U(i).name={[events{j},'3-10']};
SPM.Sess(sess).U(i).ons=event.values(3:10);
SPM.Sess(sess).U(i).dur=event.dur(3:10);
SPM.Sess(sess).U(i).P.name='r(t-1)';
SPM.Sess(sess).U(i).P.h=1;
SPM.Sess(sess).U(i).P.P=FileInfo.TruRatio(2:end-1)';
SPM.Sess(sess).U(i).P(end+1).name='i(t-1)';
SPM.Sess(sess).U(i).P(end).h=1;
SPM.Sess(sess).U(i).P(end).P=FileInfo.InvRatio(2:end-1)';
SPM.Sess(sess).U(i).P(end+1).name='r(t-2)';
SPM.Sess(sess).U(i).P(end).h=1;
SPM.Sess(sess).U(i).P(end).P=FileInfo.TruRatio(1:end-2)';
SPM.Sess(sess).U(i).P(end+1).name='i(t-2)';
SPM.Sess(sess).U(i).P(end).h=1;
SPM.Sess(sess).U(i).P(end).P=FileInfo.InvRatio(1:end-2)';
i=i+1;
else
SPM.Sess(sess).U(i).name={events{j}};
SPM.Sess(sess).U(i).ons=event.values;
SPM.Sess(sess).U(i).dur=event.dur;
SPM.Sess(sess).U(i).P(1).name='none';
i=i+1;
end
end
if length(motion) == 0, % no movement parameters
SPM.Sess(sess).C.C=[];
SPM.Sess(sess).C.name={};
elseif size(motion, 2) > 20;
SPM.Sess(sess).C.name={'first scan', 'last scan'};
SPM.Sess(sess).C.C=motion;
else % assume movement parameters
SPM.Sess(sess).C.C=motion;
SPM.Sess(sess).C.name={'X','Y','Z','x','y','z'};
end
% global normalization: OPTIONS:'Scaling'|'None'
SPM.xGX.iGXcalc='None';
% low frequency confound: high-pass cutoff (secs) [Inf=no filtering]
SPM.xX.K(sess).HParam=128;
% intrinsic autocorrelations: OPTIONS: 'none'|'AR(1) + w'
SPM.xVi.form='AR(1) + w'; %'AR(1) + w';
SPM.nscan(sess)=size(f,1);
SPM.xY.P=f;
SPM.xY.RT=2; % Time between scans
SPM=spm_fmri_spm_ui_m(SPM,'');
SPM=spm_spm_m(SPM);
succeed=true;
catch
e=lasterror;
e.message
e.identifier
for i=1:length(e.stack)
s=e.stack(i);
s.file
s.name
s.line
end
end
path(paths);
cd(pwds);
warning on;
if ~succeed&&tries>1
system('sleep 10');
analyze(runid,group,subj,tries-1);
end
function ttest6(runid,imgs)
% Perform a ttest on beta images. Creates two images as output
% corresponding to ttest's h and p parameters.
%
% Purpose: Verify our beta maps satisfy King-Casas et al 2008 results.
%
% input:
% runid - unique run identifier for all groups
% imgs - which images to use (folder names)
%
% Misha Koshelev
% March 4th, 2010
% Montague Laboratory
if ~exist('runid')
runid='600';
end
source=['/mnt/nfs/proj/trust_agnostic',filesep,runid];
target=[source,filesep,'other',filesep,'ttest'];
if ~exist(source,'dir')
keyboard
end
% if exist(target,'dir')
% system(['rm -rf ',target]);
% end
if ~exist(target,'dir')
mkdir(target);
end
warning off;
% save paths
paths=path;pwds=pwd;
addpath(genpath(pwd));
img=dir(source);
if length(img)>=2
img(1:2)=[];
end
img(find(~[img.isdir]))=[];
img(find(strcmp({img.name},'other')))=[];
if ~exist('imgs','var')
fprintf('Images:\n');
for i=1:length(img)
fprintf('%d. %s\n',i,img(i).name);
end
return;
else
suffix='';
for idx=1:length(imgs)
i=imgs(idx);
if idx>1
suffix=[suffix,'_'];
end
suffix=[suffix,img(i).name];
end
end
if length(imgs)~=1&&length(imgs)~=2
keyboard
end
M=cell(length(imgs),1);
for idx=1:length(imgs)
i=imgs(idx);
f=dir([source,filesep,img(i).name,filesep,'*.img']);
if length(f)>2
f(1:2)=[];
end
M{idx}=Inf*ones(length(f),41,48,35);
for j=1:length(f)
V=spm_vol([source,filesep,img(i).name,filesep,f(j).name]);
M{idx}(j,:,:,:)=spm_read_vols(V);
end
end
if length(imgs)==1
[h,p,ci,stats]=ttest(M{1});
else
[h,p,ci,stats]=ttest2(M{1},M{2});
end
V.fname=[target,filesep,'h_',suffix,'.img'];
spm_write_vol(V,squeeze(h));
V.fname=[target,filesep,'p_',suffix,'.img'];
spm_write_vol(V,squeeze(p));
save([target,filesep,'ttest_',suffix,'.mat'],'h','p','ci','stats');
% restore paths
path(paths);
cd(pwds);
warning on;
function analyze(runid,group,subj,tries)
% Analyze investors using SPM2
%
% input:
% runid - unique run identifier
% group - group identifier
% subj - subject
% tries - # of times to retry on failure
%
% Misha Koshelev
% February 24th, 2010
% Montague Laboratory
if ~exist('tries','var')
tries=3;
end
if regexp(runid,'3.*')
analyze3(runid,group,subj,tries);
return;
elseif strcmp(runid,'4')
analyze4(runid,group,subj,tries);
return;
elseif regexp(runid,'6.*')
analyze6(runid,group,subj,tries);
return;
elseif strcmp(runid,'7')
analyze7(runid,group,subj,tries);
return;
elseif strcmp(runid,'8')
analyze8(runid,group,subj,tries);
return;
end
warning off;
group
subj
% save paths
paths=path;pwds=pwd;
addpath(genpath(pwd));
succeed=false;
try
% load data
target=['/mnt/nfs/proj/trust_agnostic',filesep,group,filesep,subj];
cd([target,filesep,'other']);
load('FileInfo');
load('motion');
load('source');
source=[source,filesep,'GLMs',filesep,subj];
cd('..');
if exist(runid,'dir')
system(['rm -rf ',runid]);
end
mkdir(runid);
% image files
cd([source,filesep,'run1']);
f=spm_get('Files',pwd,'sw*.img');
cd([target,filesep,runid]);
spm_defaults;
global defaults;
% HRF
SPM.xBF.name='hrf'; %'hrf (with time derivative)';
SPM.xBF.order=1;
SPM.xBF.length=32.1250; % length in seconds
SPM.xBF.T=16; % number of time bins per scan
SPM.xBF.T0=1; % middle slice/timebin
SPM.xBF.UNITS='secs'; % OPTIONS: 'scans'|'secs' for onsets
SPM.xBF.Volterra=1; % OPTIONS: 1|2=order of convolution
% Trials
sess=1;i=1;
events=fieldnames(FileInfo.EventStruct{1});
for j=1:length(events)
event=FileInfo.EventStruct{1}.(events{j});
if strcmp(events{j},'ScreenInvAmount') % 4359 SHOWING MESSAGE <Make Initial Investment> AND PIC <0>
% rounds 1&2 only - no regressor
SPM.Sess(sess).U(i).name={[events{j},'1-2']};
SPM.Sess(sess).U(i).ons=event.values(1:2);
SPM.Sess(sess).U(i).dur=event.dur(1:2);
SPM.Sess(sess).U(i).P(1).name='none';
i=i+1;
% rounds 3-10 - regressed with i&r's at past 2 time points
SPM.Sess(sess).U(i).name={[events{j},'3-10']};
SPM.Sess(sess).U(i).ons=event.values(3:10);
SPM.Sess(sess).U(i).dur=event.dur(3:10);
SPM.Sess(sess).U(i).P.name='r(t-1)';
SPM.Sess(sess).U(i).P.h=1;
SPM.Sess(sess).U(i).P.P=FileInfo.TruRatio(2:end-1)';
SPM.Sess(sess).U(i).P(end+1).name='i(t-1)';
SPM.Sess(sess).U(i).P(end).h=1;
SPM.Sess(sess).U(i).P(end).P=FileInfo.InvRatio(2:end-1)';
SPM.Sess(sess).U(i).P(end+1).name='r(t-2)';
SPM.Sess(sess).U(i).P(end).h=1;
SPM.Sess(sess).U(i).P(end).P=FileInfo.TruRatio(1:end-2)';
SPM.Sess(sess).U(i).P(end+1).name='i(t-2)';
SPM.Sess(sess).U(i).P(end).h=1;
SPM.Sess(sess).U(i).P(end).P=FileInfo.InvRatio(1:end-2)';
i=i+1;
else
SPM.Sess(sess).U(i).name={events{j}};
SPM.Sess(sess).U(i).ons=event.values;
SPM.Sess(sess).U(i).dur=event.dur;
SPM.Sess(sess).U(i).P(1).name='none';
i=i+1;
end
end
if length(motion) == 0, % no movement parameters
SPM.Sess(sess).C.C=[];
SPM.Sess(sess).C.name={};
elseif size(motion, 2) > 20;
SPM.Sess(sess).C.name={'first scan', 'last scan'};
SPM.Sess(sess).C.C=motion;
else % assume movement parameters
SPM.Sess(sess).C.C=motion;
SPM.Sess(sess).C.name={'X','Y','Z','x','y','z'};
end
% global normalization: OPTIONS:'Scaling'|'None'
SPM.xGX.iGXcalc='None';
% low frequency confound: high-pass cutoff (secs) [Inf=no filtering]
SPM.xX.K(sess).HParam=128;
% intrinsic autocorrelations: OPTIONS: 'none'|'AR(1) + w'
SPM.xVi.form='AR(1) + w'; %'AR(1) + w';
SPM.nscan(sess)=size(f,1);
SPM.xY.P=f;
SPM.xY.RT=2; % Time between scans
SPM=spm_fmri_spm_ui_m(SPM,'');
SPM=spm_spm_m(SPM);
succeed=true;
catch
e=lasterror;
e.message
e.identifier
for i=1:length(e.stack)
s=e.stack(i);
s.file
s.name
s.line
end
end
path(paths);
cd(pwds);
warning on;
if ~succeed&&tries>1
system('sleep 10');
analyze(runid,group,subj,tries-1);
end
function ttest6(runid,imgs)
% Perform a ttest on beta images. Creates two images as output
% corresponding to ttest's h and p parameters.
%
% Purpose: Verify our beta maps satisfy King-Casas et al 2008 results.
%
% input:
% runid - unique run identifier for all groups
% imgs - which images to use (folder names)
%
% Misha Koshelev
% March 4th, 2010
% Montague Laboratory
if ~exist('runid')
runid='600';
end
source=['/mnt/nfs/proj/trust_agnostic',filesep,runid];
target=[source,filesep,'other',filesep,'ttest'];
if ~exist(source,'dir')
keyboard
end
% if exist(target,'dir')
% system(['rm -rf ',target]);
% end
if ~exist(target,'dir')
mkdir(target);
end
warning off;
% save paths
paths=path;pwds=pwd;
addpath(genpath(pwd));
img=dir(source);
if length(img)>=2
img(1:2)=[];
end
img(find(~[img.isdir]))=[];
img(find(strcmp({img.name},'other')))=[];
if ~exist('imgs','var')
fprintf('Images:\n');
for i=1:length(img)
fprintf('%d. %s\n',i,img(i).name);
end
return;
else
suffix='';
for idx=1:length(imgs)
i=imgs(idx);
if idx>1
suffix=[suffix,'_'];
end
suffix=[suffix,img(i).name];
end
end
if length(imgs)~=1&&length(imgs)~=2
keyboard
end
M=cell(length(imgs),1);
for idx=1:length(imgs)
i=imgs(idx);
f=dir([source,filesep,img(i).name,filesep,'*.img']);
if length(f)>2
f(1:2)=[];
end
M{idx}=Inf*ones(length(f),41,48,35);
for j=1:length(f)
V=spm_vol([source,filesep,img(i).name,filesep,f(j).name]);
M{idx}(j,:,:,:)=spm_read_vols(V);
end
end
if length(imgs)==1
[h,p,ci,stats]=ttest(M{1});
else
[h,p,ci,stats]=ttest2(M{1},M{2});
end
V.fname=[target,filesep,'h_',suffix,'.img'];
spm_write_vol(V,squeeze(h));
V.fname=[target,filesep,'p_',suffix,'.img'];
spm_write_vol(V,squeeze(p));
save([target,filesep,'ttest_',suffix,'.mat'],'h','p','ci','stats');
% restore paths
path(paths);
cd(pwds);
warning on;
|