Toolboxes¶

diagnostics¶

diagnostics.BiasPowerloss(tc, X, c, beta, df, z, pval)¶

Calculate the approximate bias and power loss due to mis-modeling This works with the Mismodeling Toolbox described by Loh et al. 2008

Usage:

[b bias pl Pc Pe] = BiasPowerloss(tc, X, c, beta, df, z, pval)

Inputs:	tc: fMRI time course X: design matrix for multiple regression c: contrast of interest beta: (mismodeled) beta value df: degrees of freedom z: p-value calculated from ResidScan pval: cut-off p-value
Outputs:	b: updated (correct) beta value bias: bias pl: power loss
References:	Loh, J. M., Lindquist, M. A., Wager, T. D. (2008). Residual Analysis for Detecting Mis-modeling in fMRI. Statistica Sinica, 18, 1421-1448.

Update design matrix using correct model

diagnostics.ResidScan(res, FWHM)¶

Calculates P(M>=t) where M is the max value of the smoothed residuals. In this implementation the residuals are smoothed using a Gaussian kernel.

Usage:

function [z sres] = ResidScan(res, FWHM)

Inputs:	res: residual time course FWHM: Full Width Half Maximum (in time units)
Outputs:	z: pvalues sres: smoothed residuals sres_ns: smoothed residuals (non standardized)

diagnostics.add_nuisance_to_SPMcfg(Xn)¶

Adds a matrix Xn to the end of covariates of interest in

xX structure in SPMcfg.mat

Usage:

function add_nuisance_to_SPMcfg(Xn)

Inputs:	oXn: should contain ALL nuisance covariates and intercepts as in output of tor_get_physio.m

This function is automatically run by tor_get_physio.m

diagnostics.batch_efficiency(dwcard)¶

Start in directory above individual model/results directories

Usage:

function batch_efficiency(dwcard)

Inputs:	dwcard: is a wildcard for directories to probe, e.g., ‘subject*’

diagnostics.batch_t_histograms(varargin)¶

Creates page(s) of t stat histograms for each subject level contrast in set of subject level analyses using image_intensity_histograms.

Usage:

batch_t_histograms([options])

Optional Inputs:
	{analysis_dirs}: run on all contrasts in directories of cell array {analysis_dirs} (DEFAULT: use all directories in working directory containing spmT_.img files) ‘o’, ‘output_directory’:* specify output directory to contain saved .png files

diagnostics.canlab_qc_metrics1(epi_names, mask_name, varargin)¶

Calculate quality control metrics on a 4-D EPI Analyze or Nifti file

Standard CANlab usage is to enter a single 4-D ravol* for one run, and the brain mask implicit_mask.img created in canlab_preproc.m

Optional Inputs:
Inputs:	epi_names: Names of (usually one) 4-D EPI file, in cell or string, full path mask_name: Name of brain mask image, string, full path IF EMPTY: Uses implict masking (better) and calculates ghost/signal outside mask
	noplot: skip plots noverbose: skip output printout to screen printfile: followed by name of file to print output to, full path noheader: suppress printing of header (var names) in output
Missing values and basic info:
	num_images: number of images missing_vox: Voxels in mask with NaN values or zero values at every time point missing_images: Images with NaN values or zero values at every voxel missing_values: NaN or zero values in valid images / voxels. Zeros could be interpreted as values of zero in analysis, causing artifacts in results if these are actually invalid values.

Missing voxels will often be ignored in analyses in most software, but Missing images/values could cause problems

Basic signal to noise:
	perc-mean-ghost: mean signal outside the mask / mean total signal mean_snr: mean Cohen’s d (signal/noise, SNR) across time (temporal SNR) within the mask. Mean divided by standard deviation across time at each voxel, averaged. Higher is better. snr_inhomogeneity: standard deviation of SNR within the mask. Lower is better. snr_inhomogeneity95: 95% confidence range for SNR within the mask. Lower is better. rms_successive_diffs: Essentially a high-pass filtered version of SNR, expressed as a fraction of the overall mean and averaged across voxels. Lower is better. rms_successive_diffs_inhomogeneity: standard deviation of the above across voxels. Lower is better.
Left-right asymmetry:
	signal_rms_asymmetry: Voxel-wise left/right root mean square asymmetry in mean signal across time, expressed as a fraction of the mean SNR. Reflects both gross inhomogeneity and noise. Lower is better. signal_hemispheric_asymmetry: Root mean square difference between left and right hemispheres, expressed as a fraction of the grand mean signal across time. Reflects gross inhomogeneity. Lower is better. snr_rms_asymmetry: Voxel-wise left/right root mean square asymmetry in SNR, expressed as a fraction of the mean SNR. Reflects both gross inhomogeneity and noise. Lower is better. snr_hemispheric_asymmetry: Root mean square difference between left and right hemispheres, expressed as a fraction of the mean SNR. Reflects gross inhomogeneity. Lower is better.
Examples:

%SETUP
mydir{1} = '/Users/tor/Documents/Tor_Documents/Coursework_and_Teaching/psyc7215/Data/UM_Face_House/060518mw/Functional/Preprocessed/run_01';
wcard = 'rarun*img';
epi_names = filenames(fullfile(mydir{1}, wcard), 'absolute');

maskdir = '/Users/tor/Documents/Tor_Documents/Coursework_and_Teaching/psyc7215/Data/UM_Face_House/060518mw';
mask = 'implicit_mask.img';
mask_name = fullfile(maskdir, maskname);

%RUN
QC = canlab_qc_metrics1(epi_names, mask_name);

QC = canlab_qc_metrics1(epi_names, mask_name, 'noplot', 'printfile', 'test_qc.txt');
QC = canlab_qc_metrics1(epi_names, mask_name, 'noplot', 'printfile', 'test_qc.txt', 'noheader');

diagnostics.check_cluster_data(cl)¶

loads the first 5 images from the first voxel

Usage:

check_cluster_data(cl)

cl(1).imnames(1:5,:)

diagnostics.compare_subjects(varargin)¶

This function compares a set of images to one another and does some diagnostics on the similarity among images. - It returns multivariate distances and dissimilarities among images - It works on the GLOBAL signal after standardizing each image (case 1) or the REGIONAL values in each cluster (case 2) - You can also enter a reference image, in which case each image will be correlated with the ref.

Usage:

function [ds, g, mystd, d, d2, c, c2, mi, b, eigv, eigval] = compare_subjects([img files or clusters], [mask], ...
                                   [plot flag], [title on figure], [standardize flag], [text labels], [ref image])

Inputs:

a list of image names to compare

OR

a clusters structure, with data to compare in timeseries field

If a mask is entered, only voxels in the mask (e.g., with value of 1) will be used. You can use this option to specify brain-only or gray-matter only voxels

textlab: optional text labels for each image, can be empty []

If a ref image is entered, each image will be correlated with the ref, and values will be saved for the correlation (plot 2 will show these values) Useful for comparing anatomical imgs with template, etc.

Outputs:

from correls with ref image are in variable “c”

ds:: multivariate distance (sim. to Mahalanobis) for each image ds is a matrix of squared distances, case numbers, and expected chi2 values (in columns in this order) rows are cases
g:: global value for each image
d:: global distance from mean image distance, or dissimilarity, is the average absolute deviation between images
d2:: matrix of distances among all images
c:: correlation between real valued voxels and mean image
c2:: correlations among all images (treating voxels as cases)
mi:: mutual information between images, with hist2.m
b:: principal component scores on correlation matrix for eigenvalues > 1
eigv:: eigenvectors
eigval:: eigenvalues

Examples:

% Compare normalized anatomcals with standard brain
P = get_filename2(['sub*\Anatomy\nscalped_ft1.img']);
[ds, g, mystd, d, d2, c, c2, mi] = compare_subjects(P, which('brain_avg152T1.img'), 1, 'intext_countloc', 1, [], which('avg152T1.img'));

diagnostics.compare_subjects256(varargin)¶

This function compares a set of images to one another and does some diagnostics on the similarity among images. - It returns multivariate distances and dissimilarities among images - It works on the GLOBAL signal after standardizing each image (case 1) or the REGIONAL values in each cluster (case 2) - You can also enter a reference image, in which case each image will be correlated with the ref.

Usage:

function [ds,g,mystd,d,d2,c,c2,mi,b,eigv,eigval] = compare_subjects256([img files or clusters],[mask], ...
                               [plot flag],[title on figure],[standardize flag],[text labels],[ref image])

Inputs:

a list of image names to compare

OR

a clusters structure, with data to compare in timeseries field

If a mask is entered, only voxels in the mask (e.g., with value of 1) will be used. You can use this option to specify brain-only or gray-matter only voxels

textlab: optional text labels for each image, can be empty []

If a ref image is entered, each image will be correlated with the ref, and values will be saved for the correlation (plot 2 will show these values) Useful for comparing anatomical imgs with template, etc.

Outputs:

from correls with ref image are in variable “c”

ds:: multivariate distance (sim. to Mahalanobis) for each image ds is a matrix of squared distances, case numbers, and expected chi2 values (in columns in this order) rows are cases
g:: global value for each image
d:: global distance from mean image distance, or dissimilarity, is the average absolute deviation between images
d2:: matrix of distances among all images
c:: correlation between real valued voxels and mean image
c2:: correlations among all images (treating voxels as cases)
mi:: mutual information between images, with hist2.m
b:: principal component scores on correlation matrix for eigenvalues > 1
eigv:: eigenvectors
eigval:: eigenvalues

Examples:

% Compare normalized anatomcals with standard brain
P = get_filename2(['sub*\Anatomy\nscalped_ft1.img']);
[ds,g,mystd,d,d2,c,c2,mi] = compare_subjects256(P,which('brain_avg152T1.img'),1,'intext_countloc',1,[],which('avg152T1.img'));

diagnostics.displayme(mm, txtlab, tlab2)¶

Used in img_hist2 - included as internal function there. this function is for indepenent re-display after img_hist2 is finished.

Usage:

function [subjM,Mtotalv] = displayme(mm,txtlab,tlab2)

Example:

% TO run:
[O.subjM,O.Mtotalv] = displayme(O.m,txtlab,'MEANS');
[O.subjS,O.Stotalv] = displayme(O.s,txtlab,'STD');
[O.subjW,O.Wtotalv] = displayme(O.w,txtlab,'SKEWNESS');
[O.subjK,O.Ktotalv] = displayme(O.k,txtlab,'KURTOSIS');

diagnostics.ellipse(x, v1, v2, c, varargin)¶

Gives x and y coordinates for an ellipse, given x coordinates, at a distance of c

Usage:

[x,y] = ellipse(x,v1,v2,c,[sorting method])

Inputs:

Based on the formula for an ellipse, x^2/v1^2 + y^2/v2^2 = c

c:: is the distance from the origin
v1:: is the x half-length
v2:: is the y half-length
x:: is a vector of points covering the x coordinates in the ellipse

Sorting methods:

sort by x, produces elliptical line in plot
sort by y, produces horizontal lines in plot

Examples:

[x,y]=ellipse((randn(1000,1)),1.5,2.5,1); figure; hh = plot(x(2:end-1),y(2:end-1),'r-')
rotate(hh,[0 90],45)  % rotate around z-axis by 45 degrees
x2 = get(hh,'XData'); y2 = get(hh,'YData'); hold on; plot(x2,y2,'bx');
rotate(hh,[0 90],-45)  % rotate original ellipse back

% fill
fill(x,y,'r','FaceAlpha',.2)

diagnostics.fft_calc(dat, TR)¶

Simple function to calculate the FFT power of a data vector (dat) as a function of frequency, given a sample-to-sample repetition time (TR)

Usage:

[myfft, freq] = fft_calc(dat, TR)

diagnostics.fmri_mask_thresh_canlab(fmri_file, outputname, implicit_masking_method, plotfigs)¶

Implicit determination of which voxels are in-brain, based on the intensities of functional images. Assumes much (most) of the image has near-zero background noise values, and the in-brain values are substantially higher.

Usage:

[mask_thresh, cl, inmaskvox, in_mask_logical_vector, maskfilename] = fmri_mask_thresh_canlab(fmri_file, outputname)

Implicit_masking_method:
Inputs:	fmri_file: is either a list of file names or an fmri_data object File names: a (preferably) 4-D file of imaging data, Analyze .img or .nii fmri_data object: With multiple images loaded with no mask outputname: is a mask file output name, e.g., ‘mask.img’, with .img extension. Empty [] means do not write output image.
	mean: take the top 95% of voxels above the mean value. used by default if no value is entered dip: smooth the histogram and take the top 95% of values above the first positive gradient **plotfigs [1/0]: enable or suppress mask display and orthviews
Outputs:	mask_thresh: signal-value above which voxels are considered in brain c1: clusters, from iimg_indx2clusters inmaskvox: number of inmask voxels dat: binary matrix of voxels that are in (1) or out (0) of mask

Note: we want to be more inclusive than not at this stage.

last edited Oct 2011 - add support for fmri_data/image_vector objects added figure suppression, SG 12/14/15 defaults

diagnostics.get_filename(dwcard, wcard, varargin)¶

Usage:

function P = get_filename(dir_wcard,file_wcard,[verbose])

Start in directory above individual subject directories

Inputs:	dwcard: Enter dwcard for wildcard of directories to look in; e.g. ‘02’ wcard:* Enter wcard for image or file to get - e.g., ‘beta_0001.img’ This can also include subdirectories (no ‘s) ‘anatomy/betaimg’
Output:	Returns list of all files in individual directories in string matrix

Missing files, or entries in directory that do not contain files, are removed from the list.

NOT entering a * in dwcard seems to produce an error.

Examples:

P = get_filename('02*','beta_0001*')
P = get_filename('02*','beta_0001.img')
P = get_filename('02*','anatomy/nnhe*_seg1.img')
P = get_filename('020515sp*','anatomy/nnhe*_seg1.img')

one * is allowed in first string, multiple *‘s in second, as long as they are in the filename, not directory names!

diagnostics.get_filename2(dwcard, varargin)¶

Usage:

function [P,P2,d] = get_filename2(search string (as with ls command),[verbose])

Start in directory above individual subject directories

Inputs:

dwcard:: Enter dwcard for wildcard of directories to look in; e.g. ‘02*’
wcard:: Enter wcard for image or file to get - e.g., ‘beta_0001.img’ This can also include subdirectories (no *‘s) ‘anatomy/beta*img’

Outputs:

Returns list of all files in individual directories in string matrix

P:: file names with full paths
P2:: file names only
d:: list of directories searched for files

Missing files, or entries in directory that do not contain files, are removed from the list. NOT entering a * in dwcard seems to produce an error.

Examples:

P = get_filename2('02*/beta_0001*')

one * is allowed in the directory structure right now, multiple *‘s in the filename.

diagnostics.hist2(A, B, res, varargin)¶

2-D histogram with res bins

Usage:

[H,mi,H2] = hist2(A,B,res,[plot])

A and B can be 3D, as in image volumes mi is mutual information, a la spm_mireg.m

diagnostics.image_intensity_histograms(fout, imgs, varargin)¶

Makes a sheet (fout.png) of intensity distribution histograms of imgs.: Will put an even number of rows of subplots on each page saved. If more than one page is needed, will title outputs fout_1.png, etc.

Usage:

image_intensity_histograms(fout, imgs, [options])

Inputs:

fout:: imfilename to be saved (‘.png’ will be appended)
imgs:: ima cell array of filenames of images to make histograms of

Optional Inputs:

obj:: im
‘titles’, cellarray:: use strings in cellarray as plot titles (DEFAULT: use file names from imgs)
‘bins’, n:: use n bins in histograms (DEFAULT: 100)
‘ymax’, y:: use y-axis from 0 to y (DEFAULT: 10,000)
‘xmax’, x:: use x-axis from -x to x (DEFAULT: 10)
‘cols’, c:: use c columns of subplots (DEFAULT: 5)
‘maxrows’, r:: use no more than r columns of subplots per page (DEFAULT: 7)
‘includezeros’:: include zero intensities in histograms (DEFAULT: exclude zeros)

Tor Wager .. add path if necessary

diagnostics.img_hist(imgname, subdir)¶

A general function for plotting histograms of any image For each subject, comparing across subjects

Inputs:	imgname: name of image file to make intensity histograms from subdir: cell array of text strings containing names of individual subject directories (wherein are contained the file specified in imgname or each subject)

Performs the histogram plot twice, once for CSF space and once for gray matter

Locations of gray and CSF masks for each subject must be defined in the defaults section of the script. (hard-coded)

Start in directory above individual subject results

Examples:

img_hist('beta_0010.img',subdir)
img_hist('con_0002.img',{'020827mk' '020829jh' '020903lb'}

% for batch
d = dir('020726ag/beta*img'); d = str2mat(d.name);
for i = 1:10:size(d,1)
    img_hist(deblank(d(i,:)),EXPT.subjects)
end

for i = 2:19, 
   if i < 10, myz = '000';, else, myz = '00';, end, 
   img_hist(['con_' myz num2str(i) '.img'],EXPT.subjects);, 
end

Tor Wager .. ..

defaults

diagnostics.img_hist2(subdir)¶

A general function for plotting histograms of any image For each subject, comparing across subjects

Inputs:	imgname: name of image file to make intensity histograms from subdir: cell array of text strings containing names of individual subject directories (wherein are contained the file specified in imgname or each subject)

Performs the histogram plot a number of times, without plotting and reports the variance in pdf moments as a function of subject, run, and condition (beta img within run).

Start in directory above individual subject results

Examples:

img_hist2(EXPT.subjects)
img_hist2({'020827mk' '020829jh' '020903lb'})

Tor Wager .. ..

defaults

diagnostics.joint_hist(x, y, varargin)¶

Create 2-D joint histogram from vectors x and y

Usage:

[z, xbins, ybins] = joint_hist(x,y,[nbins],[noplot])

Inputs:	x and y:** are vectors of paired observations on two variables
Outputs:	z: is the matrix representing the joint histogram cols of z are bins of x, rows are bins of y in plot, X axis is y, Y axis is x

Optional: number of bins, suppress plotting

Examples:

z = joint_hist(nnmfscores{i}{j}(:, 1),nnmfscores{i}{j}(:, 2), 50, 'noplot');
h = plot_joint_hist_contour(z, [0 0 1]);

diagnostics.make_conv_mtx(sz, sampres)¶

Constructs the matrix (H) for a linear convolution With the canonical SPM hrf such that Hx = conv(x,hrf)

Usage:

function H = make_conv_mtx(sz,sampres)

Inputs:	sz: size of output matrix (elements) sampres: spm_hrf sampling resolution (~ TR), OR if a vector, a custom HRF sampled at the appropriate frequency.

Tor Wager ..

diagnostics.multivar_dist(X)¶

multivariate normality checking and diagnostic plots

Usage:

[ds, S, p] = multivar_dist(X)

Input:	given matrix X with cases = rows, cols = variables
Outputs:	ds: is matrix of squared distances, case numbers, and expected chi2 values (in columns in this order) rows are cases NOTE: Sorted in order of ascending distance! S: estimated covariance matrix mv_distance: squared distances in original order of rows p: p-values in original order of rows

center

diagnostics.orthogonalize(mX, X, varargin)¶

orthogonalizes X with respect to mX, optionally scaling predictors of X For each nuisance covariate (column of X)

Usage:

function X = orthogonalize(mX,X,[scale])

Regresses out model fits and saves residuals in X

diagnostics.power_from_variance(con, N, sig2b, sig2wi, pthresh)¶

Power and effect size measures, given contrast, N, and variance component estimates

Inputs:

con:

contrast/effect magnitude estimate; “mean difference”

N:

sample size

sig2b:

between-subjects variance estimate

sig2wi:

within-subjects variance estimate: note this is not the “raw” within-subjects variance; it is the contribution to the group (2nd-level) variance, which is sig2within / number of images within-subjects

pthresh:

alpha (Type I error) rate; p-value threshold for power: calculation

con, sig2b, and sig2wi can all be vectors, so you can run this function voxel-wise for a whole map at once

Outputs:	power: Power from 0 to 1 t: effect size : expected t-value d: effect size : Cohen’s d

see effect_size_map.m for a whole-brain, image-based power mapping function

t-value threshold for significance at alpha level pthresh

diagnostics.power_loss(y, ons, X)¶: ‘true’ model fit assume ‘true’ is FIR estimate

diagnostics.publish_scn_session_spike_id(inputimgs, SUBJDATA)¶

This function is a wrapper function to call scn_session_spike_id in ‘multi-session’ mode, using input data across the runs for a single subject. It runs the program, and generates both a yaml-format text file for uploading into the CANlab database, and an html file with all the results and images for that subject embedded.

Inputs:	inputimgs: is a cell array of images (4-D) for each run in a separate cell. SUBJDATA: Input fields of SUBJDATA define the experiment name, subject name, and directories for saving both QC images + yaml and HTML
Examples:

SUBJDATA.study = 'NSF';
SUBJDATA.subject = subjects{i};
SUBJDATA.html_save_dir = fullfile(output_basedir, 'html_output');
SUBJDATA.subject_dir = fullfile(output_basedir, 'SubjectData', 'denoised_canlab', SUBJDATA.subject);

Initialize yaml file for database integration

diagnostics.qchist(images, Nbins, sparse, XLim, titles)¶

This function generates a histogram of activations from a set of statistic images. Generally, you want the images to have a normal distribution. Highly skewed distributions may be indicative of bad data.

Usage:

function: h = qchist(dat,Nbins,sparse,XLim)

  This function may generate multiple figuresa with 30 histograms
  each

Optional Inputs:
Inputs:	images: List of image file names OR fmri_data object.
	Nbins: Number of bins in each histogram (default = 100) sparse: flag for generating ONLY histograms (default = 0) XLim: Xlim (default = [-1 1]) titles: a cell array of subplot titles. If omitted, titles are inferred from assuming images come from a directory structure that looks like the following: /.../subjname/contrastimage.nii

diagnostics.reset_SPMcfg()¶: resets columns in SPMcfg by removing all non-intercept nuisance covariates. runs on the SPMcfg.mat file in the current directory

diagnostics.scale_imgs_by_csf(hP)¶

Takes a string matrix of image file names finds the mean and std of the CSF space specified in a mask (hard-coded) and standardizes images by these values

Usage:

Pout = scale_imgs_by_csf(hP)

Writes SC* images (SCaled)

assumes images are spatially normalized. uses a canonical CSF mask!

diagnostics.scn_component_rsquare(V, nuisanceX, designX)¶

Print a table of r-square values (variance explained) for each of V data vectors by nuisance (mvmt, physio) and task-related predictors

Designed to work with components

Examples:

% Typical operation
scn_component_rsquare(compscore, movement_params(1:157, :), X(1:157, :));

% No design
scn_component_rsquare(compscore, movement_params(1:157, :));

% Neither design nor nuisance, uses linear drift
scn_component_rsquare(compscore, []);

diagnostics.scn_session_spike_id(imgs, varargin)¶

Gets global image values for a session, and uses trimts.m to find outliers. The optional input MADs allows one to lower or raise the threshold for identifying scans as spikes (default = 10).

Usage:

[g, spikes, gtrim, nuisance_covs, spikesperimg, snr] = scn_session_spike_id(imgs,'mask',[mask name],'MADs',[MADs],'doplot',[0/1])

Multi-session mode returns much more output and more images, and takes in a cell array with images (preferably 4-D) for each session (run).

Inputs:	‘mask’,[pathtomaskfile]: mask images using the mask in pathtomaskfile, default: implicit mask ‘MADs’,[scalar]: change Mahalanobis distance, default: 10 ‘doplot’,[0 / 1]: plot result figures, default: true

Returns:

g:

global values

spikes:

identified spikes

gtrim:

trimmed/adjusted global values, can be used as covariate in GLM

nuisance_covs:

a matrix of 1)gtrim and 2) dummy regressors that can be used to minimize spike influence in GLM

We may want to save norms on the number of outliers found.

Examples:

% Get image names
for i = 1:6, sess_images{i} = filenames(sprintf('run%02d/vol0*img', i), 'char', 'absolute'); end

% Run
[g, spikes, gtrim, nuisance_covs, snr] = scn_session_spike_id(sess_images);

diagnostics.scn_spm_choose_hpfilter(spm_results_dir, varargin)¶

Plots and choice of optimal high-pass filter from an SPM first-level model directory (with statistics and contrasts estimated.)

Usage:

scn_spm_choose_hpfilter(spm_results_dir, ['events_only'])

SPM5 compatible and SPM8.

Called by: scn_spm_design_check.m For all regressors or events only: see scn_spm_choose_hpfilter.m

diagnostics.scn_spm_design_check(spm_results_dir, varargin)¶

Run in a single-subject (first-level) SPM directory to check design matrix variance inflation factors and high-pass filtering. Prints out table of regressors and their above-threshold VIFs (see options). Saves .png images of the key figures.

Usage:

scn_spm_design_check(spm_results_dir, varargin)

Optional Inputs:
	‘events_only’: Show plots and diagnostics for ONLY events, not nuisance covariates or other user-specified regressors. Useful when you have many nuisance covs. ‘vif_thresh’, t’: Only regressors with a VIF > t will be printed in VIF table. ‘sort_by_vif’‘: Sort regressors in VIF table by VIF (DEFAULT: order regressors as in model).

Calls: scn_spm_choose_hpfilter.m, scn_spm_get_events_of_interest.m

Examples:

scn_spm_design_check(pwd, 'events_only');

diagnostics.scn_spm_get_events_of_interest(SPM, varargin)¶

Gets events of interest.

Usage:

wh_cols = scn_spm_get_events_of_interest(SPM, varargin)

All regressors, or events only if ‘events_only’ is input as keyword ‘from_multireg’: followed by an integer, to include first n columns from the multireg R matrix as “of interest”. only works with ‘events_only’ flag, of course.

diagnostics.scnlab_norm_check(template, wanat_files, mean_func_files, subjects)¶

Compares the similarity of one or two sets of images (wanat_files, mean_func_files) to a template image and to one another (via Malanobis distance) to determine whether some images are potential outliers. This is used to check the quality of spatial warping/normalization for a group of subjects, though it could be used for other purposes as well.

Usage:

NORM_CHECK = scnlab_norm_check(template, wanat_files, mean_func_files, subjs)

Inputs:

template:: Char array with name of image of normalization template
wanat_files:: Warped (to template) anatomical file names
mean_func_files:: Names of mean functional images These images should all be in the same space/in register.
Subjs:: Optional cell array of names for each subject, for display purposes

Outputs:

A structure with metrics (NORM_CHECK)

NORM_CHECK.global_t1:

global values of first image series (wanat_files)

NORM_CHECK.std_t1:

spatial standard deviation of first image series (wanat_files)

NORM_CHECK.names_t1:

Names for columns of NORM_CHECK.norm_vs_template

NORM_CHECK.subjects:

Cell array of names for each subject

NORM_CHECK.norm_vs_template:

Similarity data for subjects (rows) x metrics (cols): {‘Dist. from group, actual chi2’, ‘Mutual info with template’, ‘Correlation with template’};

NB: Leave mean_func_files empty (e.g., []) to only check structural images

Computes metrics on the goodness of normalization based on multivariate distance, mutual information, and correlation with template. Automatically saves a .mat file of the results into the current directory.

the template file (i.e., avg152T1.nii) must be in the CURRENT working directory and have read/write permissions

USES the subfunction compare_subjects, which may be useful as a stand-alone function.

USED in canlab_preproc_norm_check.m

diagnostics.scnlab_norm_check3(wt1, subjlabels, template, mask, varargin)¶

WARNING: scnlab_norm_check3 is deprecated! All improvements are being placed in scnlab_norm_check.

Usage:

EXPT = scnlab_norm_check3(wt1,subjlabels,template,mask,[print out MI])

Inputs:	wt1: char array of wT1.img files, one per line subjlabels: cell array of subject labels template: template img that everything has been normalized to (usually the avg152T1.img file) mask: image to mask with
Optional Input:	print out mutual information table - flag for whether or not to print out the MI table; defaults to 0
Examples:

cd(studyroot); % wherever your study root is
wt1s = filenames('hr*/structural/wT1.img', 'char', 'absolute'); % assuming that hr is your study code
subjlabels = filenames('hr*');
template = which('avg152T1.nii');
mask = filenames('scalped_avg152T1_graymatter.img', 'char', 'absolute'); % set to wherever your mask is...

EXPT = scnlab_norm_check3(wt1, subjlabels, template, mask);

THIS FUNCTION IS DEPRECATED; SCNLAB_NORM_CHECK IS PREFERRED

diagnostics.scnlab_pca_check1(imgs, realign_files, X, spersess)¶

Usage:

function scnlab_pca_check1(imgs, realign_files or params (t x 6) across all runs, X, spersess)

Inputs:	imgs: list of all image names in order realign_files: movement param file for each session, names in a cell array, OR a t x 6 matrix of realignment parameters across all sessions X: design matrix; no intercept is needed
Examples:

% setup code for auditory oddball data
cd('/Users/tor/Documents/Tor_Documents/Coursework_and_Teaching/Mind_Res_Net_fMRI_Course_2008/data/auditory_oddball/2subjects-processed/s01/')

imgs = filenames('*/sw*img','absolute','char')
realign_files = filenames('*/rp*txt')

% LOAD TASK ONSETS and CREATE DESIGN MATRIX
onsets{1} = load('novel_stimuli_run1.asc');
onsets{2} = load('target_stimuli_run1.asc');
onsets{3} = load('standard_stimuli_run1.asc');
onsets{4} = load('novel_stimuli_run2.asc');
onsets{5} = load('target_stimuli_run2.asc');
onsets{6} = load('standard_stimuli_run2.asc');

regs_per_sess = 3;
nsess = 2;
for i = 1:length(onsets), onsets{i} = onsets{i}'; end
X = cell(1, nsess);
X{1} = onsets2delta(onsets(1:3), 1, 249);
X{1} = X{1}(:, 1:end-1);
X{2} = onsets2delta(onsets(4:6), 1, 249);
X{2} = X{2}(:, 1:end-1);
X = blkdiag(X{:});

diagnostics.spm_general_hist(hP, mP, textlab, varargin)¶

Usage:

function M = spm_general_hist(hP,mP,textlab,[suppress plot - enter anything])

Inputs:	hP: list of file names to compute histograms from mP: list of file names to compute masks from textlab: text string, e.g. ‘ventricles’ to label output tiffs
Output:	histograms for all input images (usually contrast images from individual subjects) plotted against a normal curve.

The expected output is that each image will have roughly mean 0, with bumps or tails in the distribution of there are real activations in some parts of the brain.

Looking at these histograms may be helpful for detecting outliers or subjects with strange contrast values. These may be caused by bad scaling, multicolinearity in the design matrix, acquisition artifacts, task-correlated head movement, or ???

Histograms (blue) are overlaid on a Gaussian distribution (red) with a mean of 0 and a standard deviation equal to that of the observed data.

diagnostics.spm_rfx_hist(cwd)¶

Usage:

function spm_rfx_hist(cwd)

Input:	a directory name where an spm or SnPM random effects analysis lives
Outputs:	histograms for all input images (usually contrast images from individual subjects) plotted against a normal curve.

The expected output is that each image will have roughly mean 0, with bumps or tails in the distribution of there are real activations in some parts of the brain.

Looking at these histograms may be helpful for detecting outliers or subjects with strange contrast values. These may be caused by bad scaling, multicolinearity in the design matrix, acquisition artifacts, task-correlated head movement, or ???

Histograms (blue) are overlaid on a Gaussian distribution (red) with a mean of 0 and a standard deviation equal to that of the observed data.

diagnostics.struct2yaml(yamlfilename, DB, yamlfilemethod, dbmethod)¶

Usage:

struct2yaml(yamlfilename, DB, yamlfilemethod, dbmethod)

Inputs:	yamlfilemethod: ‘new’ or ‘add’ (append) dbmethod: how the canlab database will handle the record. ‘add’, ‘replace’, or ‘keep_existing’

translate structure into YAML format text file this will be interpretable by the canlab database

Examples:

yamlfilename = 'YAML_tmp.yaml';

DB.study = 'NSF';     % string; study code letters
DB.subject = '001';   % string; subject ID number
DB.occasion = '21';   % string; occasion ID; unique to subj*session
DB.unique_id = [DB.study '_' DB.subject '_' DB.occasion];
DB.mean_spikes_per_image = mean(cat(2, spikesperimg{:}));

struct2yaml(yamlfilename, DB, 'add', 'replace');

diagnostics.tor_get_physio(varargin)¶

Usage:

[X,mP,spmP] = tor_get_physio([mP],[spmP],[nvoxels],[doortho])

arguments are optional, but you must enter them in this order.

Tor Wager 10/21/02

Get nuisance covariates likely to be related to physiological noise and head motion The algorithm:

The program extracts raw/preprocessed image data from the ventricles (CSF space), as defined by a mask denoting which voxels are CSF for that subject. Either all voxels or a randomly selected subset [nvoxels] is subjected to principal components analysis, to determine regular patters of drift over time and across voxels. Those patterns are expected to be related to global signal drift, head movement, and physiological noise, and are assumed to be UNrelated to the task of interest, by virtue of the fact that they occur in the ventricles.

PCA is done twice on the timeseries’ of CSF voxels. The first time, PCA is done on the sums of squared values (not the correlations) of voxel timeserieses across the entire experiment, mean-centered based on the whole experiment. Most of the coherent variation in this case is expected to be due to head movement and changes in shims/gradients/etc. from run to run. The SS values are used because we want to weight the voxels with the highest variation most heavily, as they are presumably picking up most of this signal. The first 3 eigenvariates (canonical timeseries) are saved.

Following, a separate, second PCA is done on the correlation matrix of data within each session. Session data for each voxel are mean-centered and scaled relative to the session (variance of each voxel = 1). We do this because physiological noise-related signals may produce periodic signals of different magnitudes in different voxels, and we want to extract the most coherent signals we can within each session. So these eigenvariates are expected to reflect primarily noise related to physiology (heart rate, respiration). Up to 5 eigenvariates for each session are saved (nothing with eigenvalue < 1 is saved).

Next, the CSF-related nuisance covariates (eigenvariates from PCA) are combined with existing nuisance covariates and intercept columns from the existing design matrix (SPMcfg xX). The proportion of variance in each predictor of interest explained by this nuisance basis set is calculated using regression, and the nuisance covariates are orthogonalized with respect to each predictor of interest. There are good and bad results of this step. The bad is that any signal that tracks the predictors is attributed to the task, not to noise, even if it’s actually caused by physiological artifact. So the orthogonalized basis set does not protect you from physiology or movement-related false positives. However, the nuisance covariates are also unlikely to reduce power in estimating you effects of interest. More importantly, it avoids false positives created when one predictor (A) is more highly correlated with the nuisance covariates than another (B). In practice, betas for A will tend to be smaller than B, given the same actual response to both, and a random effects analysis on A-B will produce false positive activations. Orthogonalization of the nuisance set precludes this.

Inputs:

mP:: CSF mask image file. *_seg3.img output from SPM is appropriate should be in same space and have same dims as functionals but automatic reslicing is done if necessary.
spmP:: name (full path name preferred) of SPMcfg.mat file to use This contains the design matrix and raw/preproc image file names to use.
nvoxels:: Number of CSF voxels to use in PCA analysis More than 100 can be very slow and memory intensive. Fewer than 100 voxels loads a different way, and may be slower. Best is probably between 100 - 1000. 800 runs pretty fast.
doortho:: Orthogonalize nuisance covariates with respect to regs of interest This assumes that any signal that covaries with the task is, in fact, due to the task, so it gives you some bias towards finding positive results. However, the alternative is that nuisance covariates may soak up variance related to the task, and you’ll miss activations. In addition, if some regressors are more colinear with the nuisance set, you can create false “activations” when comparing these regressors to other ones. This problem exists whether or not we choose to model nuisance covariates. One solution is to use the ortho when doing random effects analyses, as the sign and magnitude of nuisance-related activations would not be expected to be the same across subjects unless the variance was really task-related. Default is 1, or “yes, do orthogonalization.”

for functions called, see this .m file.

Examples:

% get filenames for SPMcfg files and CSF mask for each subject
cd C:\Tor_Documents\CurrentExperiments\intext2\RESULTS\model1
spmP = get_filename('sub*','SPMcfg.mat');
cd C:\Tor_Documents\CurrentExperiments\intext2\
mP = get_filename('sub*','anatomy/nscalped_f*seg3.img');
% Now run:
for i = 1:size(mP,1)  
    tor_get_physio(mP(i,:),spmP(i,:),300);  % 300 voxels
    pause(10); close all
end

Functions called:
	spm functions: spm_get, etc. timeseries2.m (for < 100 voxels) read_hdr.m (big-little endian dependent; validate for your data) timeseries3.m (for > 100 voxels; uses SPM’s image reading) reslice_imgs.m mask2voxel.m (only if ind2sub.m from Matlab is not found)

GLM_Batch_tools¶

GLM_Batch_tools.canlab_glm_getinfo(modeldir, varargin)¶

SUBJECT LEVEL input:

INFO = canlab_glm_getinfo(spm_subject_level_directory, option, [n])

Get information out of an spm subject level analysis.

Options:

(each option can be called by the listed letter or word + a number, when noted)

i’ ‘input’:

number of volumes and (first) volume name for each run

‘b’ ‘betas’ [n]:

beta names (for nth session)

‘B’ ‘taskbetas’ [n]:

beta names (that didn’t come from multiple regressors) (for nth session)

‘c’ ‘cons’ [n]:

contrast names (for nth contrast)

‘C’ ‘conw’ [n]:

beta names and weights for contrasts (for nth con)

‘v’ ‘image’ [n]:

create figure of design matrix (for nth session): (design matrix is multiplied by 100 for visibility) (works well for multiple runs)

‘V’ ‘taskimage’ [n]:

same as ‘image’, but only for task betas

‘imagesc’:

same as ‘image’, but uses imagesc: (works well for single runs)

GROUP LEVEL input:

INFO = canlab_glm_getinfo(robfit_group_level_directory, option, [n])

Get information out of a robfit group level analysis.

Options:

(each option can be called by the listed word or letter + a number, when noted)

Any of the subject level options can be used on a group level robfit analysis by prefixing ‘1i’ (output is generated based on the first input to the first robust analysis).

Ex:

canlab_glm_getinfo('second_level/model3','1iconw')

‘i’ ‘input’ [n]:: input contrasts by number and name (for nth analysis)
‘I’ ‘allinput’ [n]:: input images (for nth analysis)
‘m’ ‘model’:: weights by subject (i.e., directory containing input contrast images)
‘M’ ‘allmodels’ [n]:: weights and input images (for nth analysis)

Assumptions:

In some options, the first contrasts and group level analysis directories are assumed to represent the rest, which may not be the case.

Note:

group level options do not yet return a usable INFO struct.

GLM_Batch_tools.canlab_glm_group_levels(varargin)¶

Performs group level robust GLM analysis with robfit

sets up analysis
runs robfit
(optionally) make inverse p maps (for FSL viewing)
(optionally) estimate significant cluster sizes
publishes analysis with robfit_results_batch

Usage:

canlab_glm_group_levels([options])

Optional Inputs:
	‘s’, subjects: cell array of filenames of subject-level analysis directories IN modeldir (note: modeldir won’t be prepended to absolute paths) ‘m’, modeldir: filename of directory containing subject level analyses ‘o’, grpmodeldir: output directory name ‘c’, cov: a matrix describing group level model (do not include intercept, it is automatically included as first regressor) see help robfit note: requires specifying an output directory name note: ordering of inputs (rows) must match subjects ordering ‘n’, covname: a cell array of names for the covariates in cov ‘f’, covfile: a csv file will specify the group level model: first column, header ‘subject’, contains names of subject directories (to be found in modeldir) subsequent columns have covariates, headers are names of covariates name of covfile will be name of group analysis directory (placed in grpmodeldir) ‘mask’, maskimage: filename of mask image DSGN: will use the following fields of the DSGN structure: modeldir = DSGN.modeldir subjects = DSGN.subjects maskimage = DSGN.mask
Note:	A covfile will cause other specifications of subject, cov, and covnames to be ignored. If parameters are defined more than once (e.g., modeldir or subjects), only the last entered option will count.
Defaults:	subjects: all SPM.mat-containing directories in modeldir modeldir: pwd grpmodeldir: modeldir/one_sample_t_test cov: {} (run 1 sample t-test, see help robfit) covname: ‘groupmean’ mask: ‘brainmask.nii’
Options:	‘README’: prints canlab_glm_README, an overview of canlab_glm_{subject,group}_levels ‘overwrite’: overwrite existing output directories ‘noresults’: don’t run/publish robfit_results_batch ‘onlyresults’: just run/publish robfit_results_batch, don’t run robfit (assumes existing analyses) ‘whichcons’, [which cons] vector of contrasts to analyze (DEFAULT: aall subject level contrasts) see [which cons] in help robfit ‘invp’ [, target_space_image]: generate inverse p maps and resample to the voxel and image dimensions of target_space_image (viewable on dream with ~ruzicl/scripts/invpview) ‘nolinks’: do not make directory of named links to robust directories (using contrast names) ‘dream’: if you’re running on the dream cluster, this option will cause all analyses (e.g., lower level contrasts) to be run in parallel (submitted with matlab DCS and the Sun Grid Engine) Note: currently only works with MATLAB R2009a ‘email’, address:** send notification email to address when done running

GLM_Batch_tools.canlab_glm_group_levels_run1input(wd, c)¶: child process of canlab_glm_group_levels (see canlab_glm_README.txt for an overview)

GLM_Batch_tools.canlab_glm_maskstats(DIRS, MASK, varargin)¶

Returns a MASKSTATS structure containing data from robfitdir’s input subject level: SPM analyses.

Usage:

MASKSTATS = canlab_glm_maskstats(robfitdir, mask, [options])

output structure: MASKSTATS

COV:

subject x covariate matrix (EXPT.cov from robfit design)

COVNAME:

names of covariates in COV

MASK:

array of structs (1 per mask)

MASKFILE:

the filename of the mask used

SUB:

struct containing data from subject level images

CON:

struct contains data from contrast images

NAME:

cell array (1 cell per contrast) of contrast names

IMGFILES:

cell array (1 cell per contrast) of character arrays of contrast image filenames

(MEASURE):

subject X contrast matrix of measures (see canlab_maskstats)

COVxCON.(MEASURE):

arrays of covariate matrix X contrast means correlation results (RHO and P, see help corr())

GRP:

struct containing data from group level images

BETA:

struct array (1 struct per group level regressor) of data from beta images

NAME:

name of group level regressor

IMGFILES:

cell array (1 cell per robust directory) of beta image files

(MEASURE):

row vector (1 value per robust directory) of measures (see canlab_maskstats)

The following plots are saved in each mask’s directory in the plots directory:

contrast means by contrast (means are lines across subjects on x axis)
contrast means by subject (means are dots, lined up along the x axis by contrast)
group level betas (bar plot with group level regressors grouped by subject level contrasts)

If there’s more than one regressor in the group level model, for each regressor:

scatter plot of subject level contrast means against group level regressor

Arguments:

robfitdir:: a directory in which robfit was run contains robfit directories (e.g., robust0001) preferably contains EXPT.mat or EXPTm.mat
mask:: a filename or cell array of filenames of masks to apply to data

Options:

MEASURE OPTIONS:: (see canlab_maskstats) (DEFAULT: mean (within mask’s non-zero voxels))
‘cons’, connums:: (vector of contrast numbers) only include data from contrasts as specified by numbers
‘cons’, conname(s):: (string or cell array of strings) only include data from contrasts as specified by name
‘plots’:: make plots
‘od’, dir:: will save plots in dir (DEFAULT: robfitdir/stats_scatterplots)

GLM_Batch_tools.canlab_glm_publish(varargin)¶

Usage:

canlab_glm_publish(directory_specifications [options])

Directory Specification:
	‘s’, dirs: Generates HTML reports of the output from scn_spm_design_check for directories in cell array dires (string allowable for single dir). If a directory is a subject level analysis, the HTML will be generated for that subject in the analysis directory. Ex: canlab_glm_publish(‘s’,{‘model1/1011’ ‘model1/1014’}) If a directory contains subject level analyses, an HTML will be generated with output for each subject level analysis. Ex: canlab_glm_publish(‘s’,’model1’) ASSUMPTION: lower level analyses contain an SPM.mat file. ‘g’, dirs: For each “robfit directory” in cell array dirs, will run robust_results_batch on all contrast directories (e.g., robust0001/) (string allowable for single dir). (“robfit directories” contain robfit contrast directories (like robust0001)) EITHER directories contain EXPT.mat files (see help robfit) OR an EXPT struct is loaded in the workspace and a single directory is specified OR will do best to figure out info normally contained in EXPT.mat Ex: canlab_glm_publish(‘g’, {‘group_n35’ ‘group_anxiety_n35’ ‘group_sadness_n35’})
Note:	directory paths may be absolute or relative (to working directory)
Options:	‘t’, {[pthresh clustersize] ...}: Use the paired voxelwise_pthresh and minimum_cluster_size thresholds with which to produce robfit results maps. This option must follow immediately after a ‘g’ option (see above) and will only apply to the analyses specified in that option. ONLY applies to robfit directories (no bearing on lower level design checks) DEFAULT: {[.001 5] [.005 1] [.05 1]} Ex: canlab_glm_publish(‘g’, pwd, ‘t’, {[.001 1] [.005 10] [.05 10] [.01 25]}) ‘email’, address: send notification email to address when done running Ex: canlab_glm_publish(‘g’, pwd, ‘email’, 'ruzic@colorado.edu‘)

GLM_Batch_tools.canlab_glm_roistats(DIRS, ROI, varargin)¶: Use canlab_glm_maskstats instead

GLM_Batch_tools.canlab_glm_subject_levels(dsgnarg, varargin)¶

Performs lower level GLM analysis with SPM:

specifies model
estimates model
generates contrast images for model
creates directory with named links to spmT and con maps
publishes analyses with scn_spm_design_check

Usage:

canlab_glm_subject_levels(DSGN [options])

DSGN struct - defines the model and analysis parameters

canlab_glm_subject_levels(‘README’) to see description

Options:

‘README’:: prints canlab_glm_README, an overview of canlab_glm_{subject,group}_levels
‘dsgninfo’:: prints description of DSGN structure
‘subjects’, subject_list:: ignore DSGN.subjects, use cell array subject_list
‘overwrite’:: turn on overwriting of existing analyses (DEFAULT: skip existing)
‘onlycons’:: only run contrast job (no model specification or estimation) note: will overwrite existing contrasts note: to not run contrasts, simply do not include a contrasts field in DSGN
‘addcons’:: only run contrasts that aren’t already in SPM.mat option to canlab_spm_contrast_job
‘nodelete’:: do not delete existing contrasts (consider using addcons, above) option to canlab_spm_contrast_job
‘nolinks’:: will not make directory with named links to contrast images
‘noreview’:: will not run scn_spm_design_check
‘dream’:: if you’re running on the dream cluster, this option will cause all subjects to be run in parallel (submitted with matlab DCS and the Sun Grid Engine) Note: currently only works with MATLAB R2009a
‘email’, address:: send notification email to address when done running

Model specification and estimation done by canlab_spm_fmri_model_job

Contrasts are specified by canlab_spm_contrast_job_luka see that function for more info.

GLM_Batch_tools.canlab_glm_subject_levels_run1subject(wd, s)¶: child process of canlab_glm_subject_levels (see canlab_glm_README.txt for an overview)

HRF_Est_Toolbox2¶

HRF_Est_Toolbox2.Anneal_Logit(theta0, t, tc, Run)¶

Estimate inverse logit (IL) HRF model using Simulated Annealing Creates fitted curve - 3 logistic functions to be summed together - from parameter estimates

Usage:

[theta,HH,C,P] = Anneal_Logit(theta0,t,tc,Run)

Inputs:	Run: stick function tc: time course t: vector of time points theta0: initial value for the parameter vector

HRF_Est_Toolbox2.Det_Logit(V0, t, tc, Run)¶

Estimate inverse logit (IL) HRF model Creates fitted curve - 3 logistic functions to be summed together - from parameter estimates

Usage:

[VM, h, fit, e, param] = Det_Logit_allstim(V0,t,tc,Run)

Inputs:	Run: stick function tc: time course t: vector of time points V0: initial value for the parameter vector

HRF_Est_Toolbox2.Fit_Canonical_HRF(tc, TR, Run, T, p)¶

Fits GLM using canonical hrf (with option of using time and dispersion derivatives)’;

Usage:

function [hrf, fit, e, param, info] = Fit_Canonical_HRF(tc,TR,Runs,T,p)

Inputs:	tc: time course TR: time resolution Runs: expermental design T: length of estimated HRF p: Model type
Options:	p=1 - only canonical HRF p=2 - canonical + temporal derivative p=3 - canonical + time and dispersion derivative
Outputs:	hrf: estimated hemodynamic response function fit: estimated time course e: residual time course param: estimated amplitude, height and width info: struct containing design matrices, beta values etc

HRF_Est_Toolbox2.Fit_Logit2(tc, TR, Run, T, mode)¶

Fits FIR and smooth FIR model

Usage:

function [hrf, fit, e, param] = Fit_Logit(tc,Run,t,mode)

Inputs:	tc: time course TR: time resolution Runs: expermental design T: length of estimated HRF mode: deterministic or stochastic Options: 0 - deterministic aproach 1 - simulated annealing approach Please note that when using simulated annealing approach you may need to perform some tuning before use.
Outputs:	hrf: estimated hemodynamic response function fit: estimated time course e: residual time course param: estimated amplitude, height and width

HRF_Est_Toolbox2.Fit_sFIR(tc, TR, Run, T, mode)¶

Fits FIR and smooth FIR model

Usage:

function [hrf, fit, e, param] = Fit_sFIR(tc,TR,Runs,T,mode)

Inputs:	tc: time course TR: time resolution Runs: expermental design T: length of estimated HRF mode: FIR or smooth FIR Options: 0 - standard FIR 1 - smooth FIR
Outputs:	hrf: estimated hemodynamic response function fit: estimated time course e: residual time course param: estimated amplitude, height and width

HRF_Est_Toolbox2.Get_Logit(V, t)¶

Calculate inverse logit (IL) HRF model Creates fitted curve - 3 logistic functions to be summed together - from parameter estimates

Usage:

[h] = get_logit(V,t)

Inputs:	t: vector of time points V: parameters

HRF_Est_Toolbox2.HMHRFest(y, Runs, TR, nbasis, norder)¶

HRF estimation algorithm

Inputs:	y: Data matrix (#time points) by (#subjects) by (#voxels) Runs: Stick functions for each subject (#time points) by (#conditions) by (#subjects) TR: Time resolution nbasis: Number of b-spline basis norder: Order of b-spline basis

HRF_Est_Toolbox2.PowerLoss(modres, modfit, moddf, tc, TR, Run, alpha)¶

Estimates Power-loss due to mis-modeling.

Usage:

function [PowLoss] = PowerLoss(modres, modfit, moddf, tc, TR, Run, alpha)

Inputs:	modres: residuals modfit: model fit moddf: model degrees of freedom tc: time course TR: time resolution Runs: expermental design alpha: alpha value
Output:	PowLoss: Estimated power loss

HRF_Est_Toolbox2.ResidScan(res, FWHM)¶

Calculates P(M>=t) where M is the max value of the smoothed residuals. In this implementation the residuals are smoothed using a Gaussian kernel.

Usage:

function [p sres sres_ns] = ResidScan(res, FWHM)

Inputs:	res: residual time course FWHM: Full Width Half Maximum (in time units)
Outputs:	p: pvalues sres: smoothed residuals sres_ns: smoothed residuals (non standardized)

HRF_Est_Toolbox2.get_parameters2(hdrf, t)¶

Find model parameters

Height - h

Time to peak - p (in time units of TR seconds)

Width (at half peak) - w

Calculate Heights and Time to peak:

delta = 1/(t(2)-t(1));

HRF_Est_Toolbox2.hrf_fit_one_voxel(tc, TR, Runc, T, method, mode)¶

HRF estimation function for a single voxel;

Usage:

function [h, fit, e, param] =  hrf_fit_one_voxel(tc,TR,Runc,T,type,mode)

Implemented methods include: IL-model (Deterministic/Stochastic), FIR (Regular/Smooth), and HRF (Canonical/+ temporal/+ temporal & dispersion)

Inputs:	tc: time course TR: time resolution Runs: expermental design T: length of estimated HRF type: Model type mode: Mode
Options:	p=1 - only canonical HRF p=2 - canonical + temporal derivative p=3 - canonical + time and dispersion derivative
Outputs:	hrf: estimated hemodynamic response function fit: estimated time course e: residual time course param: estimated amplitude, height and width

HRF_Est_Toolbox2.ilogit(t)¶

Calculate the inverse logit function corresponding to the value t

Usage:

function [L] = ilogit(t)

Output:	L: exp(t)./(1+exp(t));

OptimizeDesign11¶

OptimizeDesign11.optimizeGA(GA)¶

Usage:

M = optimizeGA(GA)

outputs a pseudo-random list of condition codes that optimizes multiple fitness measures for fMRI task designs

more help can be found in ga_example_script.m and in Genetic_Algorithm_readme.rtf

OptimizeDesign11.optimizeGA_epochs(GA)¶

outputs a random-ordered list of condition #s that optimizes 3 fMRI considerations

Ways to avoid block designs: counterbalancing factor power lower limit cutoff pushes power higher
Why using avg power is better than efficiency: efficiency doesn’t account for 1/f model (altho here we just use a cutoff, the same thing) avoid transformation errors in using high pass filter efficiency is based on the sample size, determined by TR, of the model - but so is fft power...

OptimizeDesign11.optimize_rand_search(GA)¶

Just like the GA, but generates random designs each time!

outputs a random-ordered list of condition #s that optimizes 3 fMRI considerations%

Ways to avoid block designs: counterbalancing factor power lower limit cutoff pushes power higher
Why using avg power is better than efficiency: efficiency doesn’t account for 1/f model (altho here we just use a cutoff, the same thing) avoid transformation errors in using high pass filter efficiency is based on the sample size, determined by TR, of the model - but so is fft power...