Revision as of 10:39, 28 August 2008 (edit) Span (Talk | contribs) (→‎How to run IRESP) ← Older edit		Revision as of 14:20, 4 September 2008 (edit) (undo) Span (Talk | contribs) (→‎3dclust) Newer edit →
Line 137:		Line 137:
	If your analysis requires using cluster criteria, you’ll need to copy the 3dclust script into your experiment’s ttests directory. (Note: 3dclust is an afni function but we also have a convention of having a shell script called 3dclust that runs the afni function with your particular study parameters). You’ll be able to find one in most other experiment ttest directories.  You can also just look at the afni help file for 3dclust and create your own.		If your analysis requires using cluster criteria, you’ll need to copy the 3dclust script into your experiment’s ttests directory. (Note: 3dclust is an afni function but we also have a convention of having a shell script called 3dclust that runs the afni function with your particular study parameters). You’ll be able to find one in most other experiment ttest directories.  You can also just look at the afni help file for 3dclust and create your own.
		+
		+	To set up the 3dclust script, copy one from another study into your ttests directory, and open it in gedit or emacs to edit.
		+
		+	$ gedit 3dclust &
	===tableDump.py===		===tableDump.py===

---

Revision as of 14:20, 4 September 2008

Welcome to part 2 of the data analysis pathway. The lab manual will cover each step in the diagram below. While the steps under each column header must be completed in order, the column heads themselves may be done in any order or combination.

Within-group Analyses

3dTTest Script

SUMMARY: Calculates t-tests on specific model coefficients – that is to say, on your contrasts of interest. 3dttest can be used for within and between group comparisons. You’ll need two separate 3dttest scripts for this purpose.
INPUT: the z*reg+orig.* output from the reg script
You’ll also need to pull coefficient tags from the z*reg+orig file. Find the correct coefficient tag by using the 3dinfo command on the z*reg+orig output file. OUTPUT: .BRIK & .HEAD files – one pair for each contrast of interest.
RUN FROM: the scripts folder
RUN ON: everyone within your group at the same time
NOTE: Before you run the script you’ll need to create a folder called “ttests” in the experiment root directory (that is, in the expt’s broadest directory level where subject folders are located).

The 3dttest script used for within group analyses is sufficiently similar to that used for between-group analyses that we will only cover a few lines of code in this section. (To learn more about the 3dttest script, see discussion in the between-group analyses section that follows.)

Between–group Analyses

3dTtest Script

SUMMARY: Calculates t-tests on specific coefficient maps – that is to say, on your contrasts of interest. 3dttest can be used for within and between group comparisons. You’ll need two separate 3dttest scripts for this purpose.
INPUT: the z*reg+orig.* output from the reg script
OUTPUT: .BRIK & .HEAD files – one pair for each contrast of interest.
RUN FROM: the scripts directory
RUN ON: everyone from both groups
NOTE: Before you run the script you’ll need to create a folder called “ttests” in the experiment root directory (that is, in the expt’s broadest directory level where subject folders are located).

The 3dTtest script can be found in the “scripts” folder in most experiment home directories. It calculates t-tests on fit coefficients of specific regressors – that is to say, on your contrasts of interest, represented now by the output of your reg script (the z*reg+orig.* file). The 3dttest script pulls the coefficients from the reg output and compares with zero (within-group) or the other group’s data (between-group). For example, it can effectively look at whether there was greater activation at anticipation for high win (usually tagged with 1 in the lookup tables) trials versus low win trials (usually tagged with -1 in the lookup tables).

The 3dttest script outputs .BRIK & .HEAD files – one pair for each contrast. These are what you view in AFNI. They’ll be in your ttests folder a format that looks something like:

$ pvnant10c9s+tlrc.BRIK
$ pvnant10c9s+tlrc.HEAD

That is

$ contrastname_#subjs+tlrc.*

You will have a different 3dttest script for each group (within group statistics), and usually only one for between-group analyses. These scripts will call up all related model coefficients (within group or between group). On occasion for less stable models you’ll have different 3dttest scripts for different between-group contrasts (eg, for endowment and risk).

The script begins by warping each subject’s z*reg+orig file into Talairach space with adwarp. It does this by taking the set of rules you created when you warped your anatomical data, and applying the exact same warp to all of your functional data. The adwarp function can also resample voxels into a new resolution [ the –dxyz 2.0 flag, for example, resamples the input into 2mm cubed voxels]. You can edit the resample size in the 3dttest script by finding that flag and editing the number that follows.

Discussion of resampling – ask Brian.

3dttest then moves on to ttests on regressors using the command 3dttest.

Here, the 3dttest script is assessing how much of an effect each regressor has on the brain activation of interest. The t-score is derived by dividing the beta value of each regressor by its standard error term. We can then assess if these t-scores are significant or not. In other words, we can now know whether our regressors are significant predictors of the brain activation.

After completing the ttests, the 3dttest script applies the 3dmerge command to convert your t-stats to z-scores, thus creating a more intuitive stat to view in AFNI.

Before you run the 3dttest script, go to the parent experiment folder and create a ttests subfolder. The script will look for this folder and put its output files there. AFNI, in turn, looks for the folder when you select the overlay. If your files are present, you’ll be able to select your contrasts of interest!

span@dmthal riskfmri]$ mkdir ttests

Now edit your 3dTtest script appropriately.

span@dmthal riskfmri]$ cd scripts
span@dmthal riskfmri]$ emacs 3dttest

You’ll be able to pull the appropriate sub-brick # of the listing for the coefficient of your contrasts from zriskbasicreg+orig or whatever your reg script output file in the format z*reg+orig is named. You can do this by using the 3dinfo command on the file (in this case a zfile output from the risk task):

span@dmthal riskfmri]$ 3dinfo zriskbasicreg+orig.*

Here is a commented within-group 3dttest script to peruse. To see an example of a between-groups 3dttest script, go to dmthal:/data3/bipomid/scripts/3dttestantoutreg#b.

Here the bracketed number [10] after each subject is the coefficient tag we’ve pulled for that particular contrast from the z*reg+orig output of the reg script. Each contrast will have only one coefficient tag.

The –prefix flag syntax basically means that whatever entry (argument) follows it is the name of the output of that contrast’s t-test. In this case, our contrast is reward versus neutral anticipation, so our output name is trvnant3 (t for t-test, 3 for number of subjects).

-base1 is what you’re comparing all the coefficients to (here 0) -set2 is your data -session flags where you want to put your output

Here we’ve commented out (#) some contrasts we’re probably not going to use in this experiment (ant and out).

Note: be sure to edit the coefficients appropriately for each new experiment. Note that 3dttest is also run FROM the scripts directory (in the program there are lines that instruct the script to look for the subject directories).

span@dmthal scripts]$ ./3dttest

Remember to run your 3dttest scripts each time you add a new subject or a number of subjects to the group. You’ll need to ratchet up the number in the contrast file names, as well as add the subject to the list of subjects whose directories are searched by the scripts. It’s worth keeping those old files around, though: you may want to see how your data changes as you add more subjects. Keep in mind that you won’t need to re-run your reg script on everybody; only on the subject who’s been added.

ANOVA

An ANOVA tests for significant differences between means from two or more groups.

VIEWING CONTRAST DATA AND THRESHOLDING YOUR DATA

(NOTE: THIS IS APPLICABLE TO BOTH IN- AND BETWEEN-SUBJECT ANALYSES.)

VIEWING YOUR CONTRAST DATA IN AFNI.

To view your group contrasts (aka group maps) in AFNI, you’ll first need to select an individual’s brain to use as the underlay. Here we’ve chosen subject kp and we copy his axi and sag anatomical files from his folder into the ttests folder for the risk experiment:

span@dmthal riskfmri]$ cp kp012207/sag* ttests/.

Now change to the ttests directory (where you’ve placed your underlay subject’s files) in your experiment home directory and open afni:

span@dmthal riskfmri]$ cd ttests
span@dmthal ttests]$ afni

In AFNI, go to Define Overlay, and set Olay to #1 and set Ulay to #1. Set ** to 1, and # to 9. Then deselect autoRange, and enter 10.

Don’t think about why you have to set Olay and Ulay to the same thing in the Define Overlay window. There is a reason, but it is hard to explain. Briefly, using it you can threshold one dataset by another dataset.

Note that the color range bar and the slider bar that you use to select your threshold for significance are NOT coupled/related in any way, though you can go in yourself and set certain colors on the color bar to correspond to certain stat ranges.

To do so, on the color range bar note that each marking corresponds to a z-score. Right now these z scores will be rather arbitrary, but if you’ve deselected autorange and set your range to 10 (so that the brightest color equals a z of 10 and the coolest a z of -10) you will now be able to relate these bars to real z-scores (and therefore p values) of interest. As you move the slider to p values you care about (.01, .001, etc), note the z scores changing to the right. Now move your color markings so that each corresponds to a z score (for example, for a z score of 3.3 you move the marking to .33, which is simply the z-score divided by the range).

Note that ULay and Olay at the bottom RH corner of the GUI (Graphical User Interface) tell you what the z-score is at the crosshair point.

Cluster Criteria: Activation Tables

3dclust

SUMMARY: 3dclust is an afni function finds clusters of a size specified by the experimenter. 3dclust tests those clusters for statistical strength at a p value specified by the experimenter. Then outputs regions meeting both of these criteria to a .txt file. We use these .txt files to create our tables of significant activation for publications. see the next section for mor details
INPUT: 1 sub-brick (aka, contrast of interest) with z-scores, e.g., zgvnant12c, as output by the 3dttest script
OUTPUT: .txt files with a table of clusters of activation that reach your specified significance
RUN FROM: the ttests directory; other studies may have a shell script called 3dclust already in the ttest directory which you can use for reference.
RUN ON: group A versus group B output for between-group, or just group A output for within group

If your analysis requires using cluster criteria, you’ll need to copy the 3dclust script into your experiment’s ttests directory. (Note: 3dclust is an afni function but we also have a convention of having a shell script called 3dclust that runs the afni function with your particular study parameters). You’ll be able to find one in most other experiment ttest directories. You can also just look at the afni help file for 3dclust and create your own.

To set up the 3dclust script, copy one from another study into your ttests directory, and open it in gedit or emacs to edit.

 $ gedit 3dclust &

tableDump.py

This script can be found on dmthal in usr/local/bin or on the scripts page of the lab manual.

This script is for creating tables in csv format from 3dclust output.

          Example1: tableDump.py 3dclustOut.txt
          Example2: tableDump.py 3dclustOut.txt mySpecialOutputName

Input: 3dclust output from afni

Optional Input: a base filename for your output. The default is to use your cluster file name.

Output: a file called table_[yourfilename].csv which contains a csv file of the region name, x, y, z coordinates and max zScore for each cluster in the 3dclust file. You will also get a file called wai_[yourfilename].csv which contains the output from the afni whereami command.

After you get the output of the tableDump scirpt you can view it in excel.

It is VERY important that you check over the region names manually. The names simply come from the "where am i" funciton in afni which often gets things wrong. Be particularly careful about the NAcc. AFNI thinks that the nacc is much much smaller than it is and will almost always say that the NAcc activation is in the Caudate.

Note: AFNI coordinates are not always the coordinates you want to report.

AFNI is a little quirky when it comes to reporting coordinates. When you place your cursor on the brain the user interface will give you X, Y, and Z coordinates in the upper left-hand corner of the main GUI (Graphic User Interface). The numbers reported up here have flipped the signs of the x and y coordinates. That is, if the x is positive it should be negative and if negative should be positive to correctly be in Talairach coordinates. (For example, if the AFNI GUI reports x = -2, y = 3, z = 7 you would report these as Talairach coordinates x = 2, y = -3, and z = 7.) The WhereAmI function in AFNI, on the other hand, reports the coordinates in the correct Talairach coordinate space; x and y’s signs are not flipped here. However, note that 3dclust script dumps coordinate info in the flipped, AFNI format – before publishing be sure to flip the signs of your X and Y coordinates.

Conversely, if you want to take Talairach coordinates that were published in a paper and go look and see where they are in a brain, you would first have to flip the x and y coordinate signs before entering them into the AFNI goto interface. However, if you’re looking at coordinates in your raw 3dclust output from AFNI, you would not need to flip the coordinates to go look at them in AFNI as they’re still in AFNI coordinate format.

The AFNI GUI does give you some cryptic information about its coordinate system – the “RAI” designation near the coordinates tells you that the Right, Anterior, and Inferior directions are represented on the negative axes.

So in AFNI, the positive directions are LPS (Left, Posterior, Superior), but when you publish the positive directions of the Talairach coordinate system are RAS (Right, Anterior, Superior).

Cluster Criteria: Bonferonni

Brian?

Cluster Criteria: ALPHASIM

Does anybody know about this? Helps with whole-brain correction.

AlphaSim is used to do a monte carlo simulation to estimate false discovery rates.

The command we use based on the slice prescription for the gamefmri project is:

 AlphaSim -nxyz 64 64 24 -dxyz 3.75 3.75 4 -iter 100000 -pthr .001 -fwhm 4 -quiet -fast

++ AlphaSim: AFNI version=AFNI_2008_02_01_1144 (Jul 3 2008) [32-bit] ++ Authored by: B. Douglas Ward ++ default NN connectivity being used

Program: AlphaSim Author: B. Douglas Ward Initial Release: 18 June 1997 Latest Revision: 10 Jan 2008

Data set dimensions: nx = 64 ny = 64 nz = 24 (voxels) dx = 3.75 dy = 3.75 dz = 4.00 (mm)

Gaussian filter widths: sigmax = 1.70 FWHMx = 4.00 sigmay = 1.70 FWHMy = 4.00 sigmaz = 1.70 FWHMz = 4.00

Cluster connection = Nearest Neighbor Threshold probability: pthr = 1.000000e-03

Number of Monte Carlo iterations = 100000

 Cl Size   Frequency    CumuProp     p/Voxel   Max Freq       Alpha
     1       8754479    0.942821  0.00100385        482    1.000000
     2        484826    0.995034  0.00011330      62539    0.995180
     3         40988    0.999449  0.00001466      31987    0.369790
     4          4503    0.999934  0.00000215       4377    0.049920
     5           543    0.999992  0.00000032        542    0.006150
     6            60    0.999999  0.00000005         60    0.000730
     7            10    1.000000  0.00000001         10    0.000130
     8             2    1.000000  0.00000000          2    0.000030
     9             0    1.000000  0.00000000          0    0.000010
    10             1    1.000000  0.00000000          1    0.000010

AlphaSim Manual

Simultaneous Inference for AFNI Data

B. Douglas Ward

Biophysics Research Institute, Medical College of Wisconsin

File:AlphaSim Ward.pdf

Download PDF

ROI (Region of Interest) Analyses and Dumping Timecourses

The reordertcdump Script

SUMMARY: Reorders and averages each trial type time course so that you have an averaged time course for each type of trial.
INPUT: motion files (3dmotion.1D), master vectors (*.vec), functional data (*.BRIK), ROI masks (Region of Interest)
OUTPUT: .tc files (tc files are vectors with time course data in the form of normalized % signal change for each ROI)
RUN FROM: scripts directory – this is where your masks should go, as well
RUN ON: All your subjects. Subjects can be run through the script separately, as well, so there’s no need to re-run everybody if you add people.

Timecourses are sets of data that record normalized percent signal change for every voxel in a region of interest (the course, over time, of that voxel, hence timecourse). You can output curves for that voxel showing the voxel’s average activity over the course of each type of trial. That is, it will pull all anticipation periods for a certain type of trial, average those together, and spit out an averaged curve for anticipation on that type of trial. It then does the same thing for the next TR (for example, choice presentation), spits out an averaged curve for choice presentation on that type of trial, and glues this on to the end of the previous, eventually giving you your averaged curve for that voxel over the course of an entire single trial type.

To “dump” timecourse data, you will first need to borrow or create masks that tag ROIs (Regions of Interest). Examples of ROIs commonly used in our lab are the mpfc, nacc, and vta. These files will be in a form that looks something like “nacc8mm+tlrc”.

If you have already created masks for another study (eg, nacc8mm – a good bet) that you would like to use in your current study, you can copy them into your new experiment’s scripts directory.

If you have no masks, you’ll need to make new ones.

I. Making a Mask

• In the scripts directory of the experiment, copy in sag116+* files, and also may need to copy in unsamp mask (mask that hasn’t been resampled to 2mm)
• Also copy in an old mask if you want to write over it
• Open AFNI. Ensure that the 'View ULay Data Brick' is selected in the Define Datamode menu. Also make sure that you’re viewing the data in Talairach view.

II. Now, in AFNI → Define Datamode → Plugins → Drawdataset (which opens an 'AFNI Editor' window), then do the following:

Copy Dataset, then Choose Dataset. It doesn't matter what old mask you copy/overwrite. (If there are no available datasets 
to copy, find an old ‘unresampled’ (normal 4mm-ish res) mask in another study’s scripts directory- they will usually end 
in xxunsamp+tlrc.)

Set these values: 
Mode: 3D Sphere R = Radius; 4mm is normal

Value: Starting from 1 for the first ROI  

Label: Usually 4-5 characters

Find the center of the region where you want to place the sphere (it is possible to zoom in and pan to get correct placement by 
using the buttons on the sides of the AFNI sagittal / coronal / axial views).  When you are centered, shift + click on the spot
to place the sphere.  The overlay that you have copied now has a region of voxels with the value as defined in the draw dataset 
box - in our case, usually a sphere with lots of 1's.

Set the Value to 2, change the Label, and repeat these steps until you’ve created all the ROIs you want.  
Note that L and R nacc  et cetera will need separate spheres that are numbered differently.

When done, click the Save As button and name the mask. 

The mask just created must be re-sampled to match the resolution of the functional data instead of the anatomical, so run 3dresample at the command line as follows:

$ 3dresample -dxyz DX DY DZ -prefix [output resampled mask name] -inset [input mask]

e.g.

$ 3dresample -dxyz 2.0 2.0 2.0 -prefix maskname -inset maskunsampfile+tlrc

Maskname is your final output mask to use in your timecourse dump.

Note that the syntax here is a bit different. Adwarp’s dxyz flag takes only one argument, here you have to specify each direction.

note: If you open an old mask in order to add new ROIs, it won't usually work (unless you do unknown tricks). So, you have to make a whole new mask with the old ROIs and the new ones.

Back to the reordertcdump script.

Now open the reordertcdump script, which can be found in your scripts directory. (If it’s not there, copy it from someone else’s directory.)

The script utilizes as input the motion files, master vectors, functional data, and the ROI (Region of Interest) masks.

Reordertcdump takes the functional data and averages the trials as described above, created one single average trial percent signal change for every single voxel, viewable as a list of numbers or as a graph. Its sole output are .tc files – columns of numbers corresponding to normalized percent signal change. Edit the script appropriately; for more detailed notes on how to do this, see the bipomid reordertcdump script located on dmthal:/data3/bipomid/scripts/reordertcdumpC.

Reordertcdump is generally run on the entire group of subjects, however, if you add a new subject later you can go back and run the script on that subject alone.

The script 'reordertcdump' (which can be copied from another exp.'s scripts dir) needs to be edited for each particular study and mask. It takes the +orig functional files for each subject, warps to Talairach space, and then extracts the time-course data from each region in the mask. Instructions on how to customize the script are given in the reordertcdump comments present in most tasks’ reordertcdump scripts.

The script first takes each subject's master vector and reorders the event types from the run. It is important to set the correct number of 'trailing TRs' so that the lagging BOLD signal resulting from each trial can be captured and plotted; currently this is set to 10 TRs.

In the script, the important function is 3dmaskave, which has the following syntax:

$ 3dmaskave -mask ${maskfile}+tlrc -mrange n n ${reordered}+tlrc > {subject}label.tc

Make sure the labels for each numbered ROI match the region masked (though not necessarily the label used in drawing the dataset in AFNI), and make sure to remember these labels for use in arrangetcs. '-mrange n n' takes values from the voxels drawn in the mask from the just-created reordered file as 'n's, averages all the voxels in the region together to produce one number for each time point, and returns a .tc output file which displays all the time points for each event type sequentially.

Note that if you need to check what values are set for a specific region, you can open the mask (open AFNI in the scripts directory) and set the overlay to the mask file. Go to the region and see what value is written there (look at the AFNI display where z-scores are normally displayed).

VIEWING TIMECOURSES AS ULAY.

It’s sometimes useful to view your group-averaged timecourses as ulay. To do this, set groupreordertcs+tlrc –type file as Ulay.

For more information on how to do this, see the commented reordertcdump script.

ALTERNATIVE METHOD OF TIMECOURSE DUMPING: IRESP

About IRESP

[Lab members, please edit as necessary; also, some of these examples are found in the endowment study; this may need to be updated as more and more newly collected studies become relevant.]

Another method of dumping out trial-type averaged timecourses was developed to deal with random ITIs. The lab’s original timecourse dumping method, using the python script reorder_regress_lag.py, was designed to dump out data from the MID task, which has a fixed – but jittered near 2 sec – trial length (*note that this also creates problems for the generation of subject mastervectors via txt2master.py, but it can handle – with a few modifications – random 2, 4, and 6 sec ITIs). Recent tasks have very random ITIs (futself) or random within a range (deldis, facetrait, endow, simplerisk). We STRONGLY recommend using pseudorandom ITIs, so to allow looking at timecourses from these tasks, a new method has been developed.

This takes advantage of AFNI’s built in impulse response (iresp) estimation capacity in 3dDeconvolve. 3dDeconvolve is used to estimate the fit of brain activity to an idealized model, but it also has the capacity to estimate iresps as well, if given a list of even onsets instead of the usual HRF-convolved regressor files.

Usage (please copy in actual script from /datae2/endowment/scripts/irespdiff):

3dDeconvolve -…….. -stim_file 1 3dmotion12.1D[1]…. …. -stim_file 7 trialtype1onsets.1D -iresp 7 -stim_maxlag 7 10 -nofout

rm –f Decon+* rm –f Decon*x1D

First, motion files are necessary so that AFNI can regresss out the effects of motion on the iresp estimation. stim_file 7 simply tags the onsets of trialtype1 with 1s and is 0 otherwise. The iresp flag tells AFNI that this is a stimfile to use for impulse response estimation and not for the normal whole-brain regression. Finally the stim_maxlag flag tells AFNI how long we would like to average activity after the event onset. A value of 10 will give a timecourse file that is 11 TRs long, counting the onset TR.

The output (besides a useless whole-brain regression output, Decon+orig) will be an AFNI 3d+time dataset containing the average response to trialtype1, separately computed for each voxel in the brain. This file can be viewed in the AFNI GUI (Graphic User Interface) and graphed if it is set to the underlay.

This file is also input into 3dmaskave (see /data2/endowment/scripts/3dmaskave), which is equivalent to the latter half of the old reordertcdump script, where the input is the trial-averaged 3d+time dataset and a mask, and the output are VOI (volume of interest) averaged timecourses as a column of 11 numbers.

Currently these iresp scripts take the output of numerous individual 3dDecons and paste them together so that all trialtypes can be dumped out (using 3dmaskave script) at one time. This is accomplished very simply by using 3dTcat (see the end of iresp scripts for an example). Then the individual trial-averaged 3d+time datasets are removed because they take up a lot of space (*this is important*).

As a way to see timecourses of particular events all over the brain, instead of using masks as spotlights on activity, advanced users may wish to take each individual’s iresp output and create a group average. This can then be plotted underneath group statistical maps in the ttests directory. This allows a check of what all areas of the brain are doing when presented with trial types of varying difficulty, for example, and can be very illuminating. These files can be created using the tcaverage script in that study. In the endowment ttests directory, various ways of dividing the data into trialtypes exist (difficulty, percent retail price, purchase decision), and the files are often broken up into trial categories (e.g. difficulty timecourses for buy trials, starting with low, then medium, then high difficulty, all in one 3d+time dataset). In this directory there is also a file that represents the average brain response to *all* trial types, representing the brain response to making a decision about a product.

A draft technical report comparing iresp timecourse dumping to reorder_regress_lay timecourse dumping is saved on the macbook insula. The first comparison looks at the response evoked by a flashing checkerboard in younger and older subjects in the ageMID dataset. These timecourses were dumped out from single voxels (*note: you can save a graphed timecourse in the AFNI GUI using options on the lower right of the graph window*) and a comparison reveals no differences. This is a simple case, but ideally both give the same results - both methods covary out motion, and this is a task with only one trial type. A comparison using a more complicated task with many trialtypes was started, using the futself study data. Here, it appeared that using iresp with many input trialtypes inside the same 3dDeconvolve program gave different – weaker – results than the original reorder_regress_lag method. Iresp may produce weaker results when it is covarying out the effects of other trial types (however, in all our designs, there should be no dependencies of trial types on one another, making the need to covary null), or it may be using an inaccurate method of estimation. Thus to get the best estimate of timecourses using iresp, we now only input one trial type onset vector into each call to 3dDeconvolve.

(*Minor note: if someone can figure out the new AFNI flag for not giving an output Decon+orig file, that would be awesome, otherwise these iresp 3dDecon calls waste time while afni computes the fits of the motion regressors, etc; also throw in the no x-matrix option, too.)

(dy 10/10/07 I think using the –nobucket flag will remove those pesky decon files)

How to run IRESP

To dump timecourses using iresp, you'll first need to make sure you have the "starter" files your scripts will need:

1. m2v lookup tables

These are not the same as our traditional m2v lookup files, which mark EVENTS of interest WITHIN trials. These mark the beginning of trial types, instead, which iresp uses to know where to start dumping those trial types for the VOI dumps. Place these m2v files in a directory called "iresp" in your expt home directory. You'll also place your iresp script here. For more on the iresp script, see below.

2. runs.1D file -- located in a folder called "vectors," the first line of this file should be 0, the second line should be the total # of TRs (end index)

3. sag116 file

4. normed, filtered data, as output by your process script (eg, norm12f+orig.BRIK)

1) Run the irespbipomid script to

• Create individualized vector .1D files using the m2v lookup tables.
• Use 3ddeconvolve to remove the contribution of motion to your data, and to create, for each subject, an average response to each specified trial type. "Impulse response" refers to the fact you flagged the _beginning_ of each trial type. It then dumps the response from that point in time forward as you specify (in this case, 10 + 1 TRs). 3dDeconvolve takes as input the .1D files created earlier in the script, as well as your runs.1D file (which specifies block length) and your normed, filtered data.

- 3dDeconvolve will output BRIK files for each trial type of interest for each individual subject. - Use 3dTCat to concatenate all the BRIK files into a single BRIK with all the trial types concatenated one after the other in the order specified. "-prefix" will specify the name of this new, uber-BRIK.

• This script can be run on individuals alone as you add them.
  • The irespbipomid script usually lives in the iresp directory with the iresp m2v files.
    • The irespbipomid script is doing FULL BRAIN dumps -- there's no ROI timecourse data yet. These are volumes with information about the average response in each and every individual voxel in the brain for that individual.

2) Then run the masterTCdump script. This script actually pulls the timecourse data from the voxels in the ROIs and averages it using 3dMaskAve. It should be run on all your subjects within the group of interest.

(a) This script takes the following files as input: your sag116 file, your mask files, and your irespbipomid BRIK file, as created by the irespbipomid script, above.

3) Now, masterTCdump drops the timecourse files into a timecourse folder in your home directory.

(a) Move one of the MatlabDemo.m files into the timecourse folder. You can find it in the old bipomid timecourse folder.

(b) Open matlab (type matlab at the prompt while in the folder), double-click on the matlabdemo file on the LH side.

(c) Open the matlabDemo.m file if it doesn't automatically open in an extra window.

(d) The matlabDemo file uses tcpaste to concatenate your subjects into groups, eg CTR & BPD, sorted by region, and concatenated in the order you specified in masterTCdump. You shouldn't run the entire file, however, as it contains different script snippets. All you need do is highlight three lines: subjects = , regions = , and tcPaste(subjects, regions). Once highlighted, right click and hit "Evaluate Selection."

4) You should now be able to open your concatenated data by double clicking on, say, "nacc_allTRctrs.csv" and importing it into your workspace. Then graph to your heart's desire.

PLOTTING TIME COURSES IN EXCEL

On insula or the windows laptop, open a template file for the given study or steal a template TC excel file from a similar study - in data/studyname/timecourses - so that graphs appear automatically. Or make a new one. Then open a .csv file for an ROI, copy all the subject's TCs and averages/stddev to the template or new Excel file, check that the graphs are plotting from the right data (select a graph, then go to graph→choose dataset) and see what you've got. It’s useful to plot the averages and stddev on one graph or sheet in the excel file, and plot each individual on the same page as the data is copied onto, so that an individual can be matched to a timecourse.

note2: Because we now use a different script to generate the baseline files in the initial processing, we no longer need to set a baseline file in reordertcdump or divide the TC data by a baseline in arrangetcs, where we now use the python script 'tcourse2excelNoBaseDiv.py' instead of the older one, 'tcourse2excel.py.'

note3: see Jamil's sweet pictures from the 'processing.presentation' in Spanlab docs for more on reordering events.

Exploratory Excel timecourse plotting is an essential part of any study. Using this, we determine what direction brain activation is going in different trials and timepoints in a priori regions of interest or those that appear in the group maps. This grants much better understanding of the data than maps alone provide. For example, if depressed subjects show greater activation than controls in the anterior cingulate when anticipating high vs. low gain rewards, what is actually going on in the ACC? Timecourses show that instead of higher activation for high gain trials in depressed subjects, the difference is due to much lower low-gain ACC activity for the depressed subjects. Thus the difference between hi and low gain trials between groups could be due to at least two underlying causes (a third would be lo>hi gain activity in the control group ACC) that are only established by plotting actual activation.

SPSS.

After you’ve dumped your timecourses to Excel .csv files, you can open them in SPSS and perform additional statistical tests (e.g., interactions of valence, differences in time courses of high versus low gain or any other contrast) that may or may not integrate behavioral information, as well.

DELTAGRAPH.

Delta Graph is a program that does graphs. It has a variety of spreadsheet and graph files. Must click on a graph to use anything in the graph file menu.

When you have a delta graph graph, and you want to copy it to somewhere else, in one of the menus you can group the graph into one item before copying.

SANITY-CHECK: MAPS AND TIMECOURSES

One of the most important things to do with freshly processed data is an integrity-check. There are several avenues you want to explore:

MAPS ADD MORE INFO ON CHECKING MAPS LATER

Red flags when looking at the output of within and between group 3dttests include:

- Unusual patterns of activation in a task that’s well-known to produce consistent, specific patterns of brain response. - Maps with enormous, uniform clusters of activation in only one direction across the entire brain.

TIMECOURSE DATA

It’s a good idea to plot your timecourses and standard errors from your Excel .csv files as soon as you’ve dumped them from the arrangetcs script. This allows you to check for outliers – individuals who may have moved during the scan, for example, or who have an incorrectly set gain or some other data processing error (this will show up as a large standard error and likely an unusual looking graph even at the group level). The latter type of error is most likely to be systematically present in the data, so unless your subjects were processed with different process scripts you can probably assume that a single subject outlier is due to that subject’s movement or a scanning error for that scan.

ADDING NEW SUBJECTS TO YOUR ANALYSIS

When adding a new subject to an existing analysis, you needn’t rerun every script on every subject. Here’s what you will need to do:

1. Construct that subject’s anatomicals with to3d.
2. Process that subject’s functional data with the process script.
3. Nudge & Warp your subject’s brain.
4. Use the txt2master script to create a .vec and .outkey file for that subject.
5. Run the reg script – it can be run on that subject only – no need to rerun for everybody.
6. Ratchet up the number of people in all your 3dttest scripts and save new versions of the scripts with those numbers in the script names. Edit the scripts to include the new subject(s) – don’t forget to increment up script internal elements that reference the total number of subjects, as well, including the names of your outputs, which if left unchanged will overwrite your old outputs and will have the incorrect number of subjects implied by the file name. Rerun the new versions of all 3dttest scripts.
7. Rerun 3dclust on the output of the 3dttest scripts.
8. There’s no need to rerun reordertcdump on everybody – you can just run it on your new subject(s) alone.
9. The arrangetcs script, on the other hand, must be run on everybody since it’s averaging timecourses across subject.

ADDING NEW REGRESSORS TO YOUR ANALYSIS

When adding a new regressor you need to rerun everybody on every script from reg onward. (Except that you wouldn’t need to do the timecourses again.)

ANALYZING NEW GROUPS IN YOUR ANALYSIS

Let’s say you’re doing a study on a patient population and controls. You may want to analyze patients on meds versus those not on meds, which will require creating new contrast groups for your analysis. When creating new contrast groups you’ll need to rerun everything from the 3dttest scripts onward, including time courses if doing group timecourse dumps.

LOGISTIC REGRESSION ANALYSIS

Stata notes

Stata is a statistical analysis program that runs on Stanford servers (e.g. Elaine). The lab primarily uses Stata to run logistic regressions, using neural activity to predict binary behavior like purchasing decisions or gamble choices. The interface is graphical, but has a convenient command-line option and easy to specify fixed-effects that make it better than SPSS. Stata is not free, however, and has not been installed on SPANlab servers. Thus lab users transfer data files to their space on a server and run Stata by ssh connection [*note: Stata should probably be pronounced like the beginning of sta-tistics]

File types: easiest to use are comma-separated (.csv) files, because Excel readily saves in this mode (ignore it’s annoying warnings) and more importantly because it is also easy to turn timecourse data into .csv files used in logits (via $ paste –d , column1.tc column2.tc.. > output.csv). Note that once a .csv file has been imported for use, Stata requires a restart to accept another import. Stata natively uses .dta files, and imported files can saved in this format for easier loading.

File structure: for logistic regressions, typical spreadsheets used in the lab have a row for each decision and columns for the various behavioral and neural variables taken from that trial (e.g. nacc activity during product presentation, product preference, price, reaction time, and choice, during trials 32 / for trial 32. Files also will have a column for subject number (best to use scanner number) and choice number.

logit Computes a logistic regression, using predictor variables(s) to estimate a binary outcome variable, e.g. logit purchase nacc will predict purchase decisions using nacc activity.

Fairly decent help can be found for this command (and others) here [1]

logistic Computes a logistic regression similar to the logit command, but only gives odds-ratios instead of t/z-scores.

regress Computes a traditional regression, predicting a scalar quantity. An example would be using preference ratings to predict product-evoked nacc activity.

pwcorr Computes a simple pair-wise correlation between two input variables. Using “pwcorr x y, sig” will additionally give you the significance of the correlation.

xi: .. i.subj Controls for fixed effects. Using “xi:” and “i.subj” (where “subj” is a column that numerically codes subjects) as bookends to any command (like regress and logit) will control for subject fixed-effects. This has been used in all SPANlab logit analyses to date.

e.g. $ xi: logit purchase nacc mpfc insula i.subj

estat ic Computes Akaike and Bayesian information criterion (AIC, BIC) for the model just run. These are measures of model strength that are better than using simple pseudo-r2 because they look at fit and additionally take a hit from extra predictor variables.

predict res, r Computes the residuals of the model just run. Initially used in the lab for covarying out motion and using residualized brain activation in logits. Puts the residuals in a new column in the input dataset. Note “res” is just the output column heading; it could be named anything.

Stata can do a lot more, but we have had no need to learn how to do ANOVAs there, for example, because SPSS works fine for small tasks on small datasets.

Brief note on mediation:

“In psychology, a mediation model is one that seeks to identify and explicate the mechanism that underlies an observed relationship between an independent variable and a dependent variable via the inclusion of a third explanatory variable, known as a mediator variable. Rather than hypothesizing a direct causal relationship between the independent variable and the dependent variable, a mediational model hypothesizes that the independent variable causes the mediator variable, which in turn causes the dependent variable.” (Wikipedia)

Mediation effects are simple to determine when you are already running logistic regressions. If the outcome variable is binary, logistic regressions can be used for the paths between the predictor and mediator and the outcome, while the “regress” regression command can be used from predictor to mediator, if the mediator is scalar. Then simply plug in the resulting coefficients and standard errors for paths a and b above into a mediation calculation.

These pages can be used to compute mediation significance: [2] or [3]

Another useful resource, with SPSS macros is here [4].

Brief note on mediation:

“In psychology, a mediation model is one that seeks to identify and explicate the mechanism that underlies an observed relationship between an independent variable and a dependent variable via the inclusion of a third explanatory variable, known as a mediator variable. Rather than hypothesizing a direct causal relationship between the independent variable and the dependent variable, a mediational model hypothesizes that the independent variable causes the mediator variable, which in turn causes the dependent variable.” (Wikipedia)

Mediation effects are simple to determine when you are already running logistic regressions. If the outcome variable is binary, logistic regressions can be used for the paths between the predictor and mediator and the outcome, while the “regress” regression command can be used from predictor to mediator, if the mediator is scalar. Then simply plug in the resulting coefficients and standard errors for paths a and b above into a mediation calculation.

This page can be used to compute mediation significance:

http://www.psych.ku.edu/preacher/sobel/sobel.htm

or here:

http://www.danielsoper.com/statcalc/calc31.aspx

Another useful resource, with SPSS macros:

http://www.comm.ohio-state.edu/ahayes/SPSS%20programs/indirect.htm