The topic of this page can be boiled down to a simple thought question: Why can't you simply shuffle voxels as a "control" for the spatial patterns that you see in your actual results?
This is a deep, far-reaching issue that is affected by many issues in (f)MRI.
What spatial structures actually are we studying / searching for?
You can think of this question as referring to the true biological structures that we are attempting to measure via MRI (or fMRI). Below are some examples.
Gray matter vs. white matter. This is a real anatomical 'structure' and there is a fairly sharp transition between gray and white matter. So, how well can we actually measure this transition?
Topography (retinotopy, tonotopy, somatotopy). Continuous/smooth progression. Clearly, this is a "low spatial frequency" structure. So, in a sense, if you are trying to measure this type of structure, the lower frequency your results, the "better".
ROIs. E.g. functionally localized ROIs. These are more questionable in terms of, are they real biological structures. But, it is still an interesting exercise to think about how well we can measure their spatial structure.
Cortical layers.
Cortical columns.
Of course, there could be many scales of spatial structures that could exist simultaneously in the brain.
MRI physics / MRI acquisition
In terms of how MRI data are acquired, there are a variety of considerations to take into account that influence the actual spatial resolution of the data you get.
Bottom line: In an ideal geometric sense, voxels are cubes. But in terms of actual signal acquisition, voxels are not cubes at all, but are more like blurry sampling elements (e.g. like a Gaussian).
T2* decay - In EPI, there is T2* signal decay over the course of your k-space measurement. This is imposing some unavoidable/intrinsic blurriness to your image in the phase-encode direction.
slice profile - Imperfections in slice excitation profile mean that you also have some intrinsic blurriness in the slice dimension. Ideal slices are perfect rectangles, but real slices are not quite like that. (Note: compared to other effects, this is probably a very minor one.)
Partial Fourier - This refers to omitting a part of k-space (typically, a block of the highest spatial frequencies). This has the effect of essentially blurring your image in the phase-encoded dimension.
k-space - Typically, the MRI pulse sequence is measuring Fourier components of your slice. Hence, voxels are really sinc functions (not cubes).
Physicists like to image phantoms. This makes the issues easier than real brains, since you don't have any head motion or physiological noise. By imaging phantoms, you can discover and see clearly certain effects such as intrinsic blur in the phase-encode dimension.
Practically speaking, what can we do?
You can accept the fact that your spatial resolution is only nominal.
You can control the slice dimension.
You know your phase-encode dimension, so you can think about whether the extra blurriness in that dimension is a problem.
Things beyond MRI acquisition
In real data, such as fMRI data acquired on a live human, there are additional issues that have to be taken into account.
There is intrinsic biological motion potentially (e.g. blood pulsation), which is yet another reason you can't actually get a perfectly sharp image.
There is also head motion over time, which means that we have to talk about what motion correction does and whether it fully compensates for displacements in head position. The short answer is it does not fully compensate, and that you get additional blurring due to the interpolation in motion correction.
Noise
Noise comes lots of flavors and variety! It is fundamental, because it will affect your analysis result, and specifically the spatial characteristics of your analysis results.
Some types of noise are structured (e.g. the effects of respiration on BOLD signals), other types of noise are not (e.g. thermal noise).
We also have to think about the signal (roughly defined as the thing you are trying to measure). The true signal you are looking for is often (but not always) low spatial frequency biased.
Noise varies across conditions, subjects, sessions. Hence, you have to be very careful and consider whether the effect you are finding (e.g. different spatial characteristics) is not inadvertently caused by some variation in noise.
Preprocessing
fMRI data are typically pre-processed, and you have to think about what these pre-processing steps are doing to the spatial resolution/accuracy of your data.
Motion correction - If the head moves, you have to "put it back" in order to do almost any reasonable analysis of the data. But interpolation is not magic. It's the best you can do, but the resulting image still has imperfections, and will have a bit of extra blurriness.
Group atlas space - If you use a group space, you should consider if anatomical variability across subjects is contributing to your analysis result.
Spatial smoothing - here, for various reasons, we deliberately induce extra intrinsic smoothness into the data
Upsampling - OBVIOUSLY, after upsampling, voxels are by no means independent.
Practical analysis ideas:
Extrapolation method. If you have multiple trials, we can try to extrapolate and see what data characteristics would have been obtained in the limit of infinite number of trials. See Figure 12 in this paper.
You could try your analysis on some empirical control data. For example, you could take some resting-state data. The advantage is that you are doing many exact matches: the same human, imaged using the exact same pulse sequence, etc.
One problem is that resting-state data has a lot of interesting spatial structure, and maybe this is too strong of a control. So, we might wonder is it possible to have more of a white noise human control, somehow?
You could deliberately use an incorrect design matrix in conjunction with your actual data.
You could make some idealize simulated data. This is very tricky, as you have to make strong assumptions about the nature of what type of noise to generate.
Relative comparisons - One approach is to just accept that you can't make absolute claims about the spatial characteristics of your data, but at least you can make relative comparisons (e.g. compare one brain region to another).