Hi, Christine-
You can:
A) first, concatenate all the 3D volumes of your subjects across the "time" axis, and then
B) compute the voxelwise stat (mean, stdev, etc.) for each voxel.
To do "A":
3dTcat -prefix ALL_SUBJ.nii.gz DSET_000 DSET_001 DSET_002 ...
where DSET_* are your FA-like volumes for each subject; I'm assuming each one is just a 3D volume; if that is *not* the case, you could use subbrick selectors, like: DSET_000'[0]' DSET_001'[0]' DSET_002'[0]', there instead.
To do "B":
3dTstat -mean -prefix STAT_mean.nii.gz ALL_SUBJ.nii.gz
3dTstat -stdev -prefix STAT_stdev.nii.gz ALL_SUBJ.nii.gz
... and see the full list of stats that can be computed in 3dTstat.
--pt