Dear Will,
thank you very much for the detailed information!
Best, Carsten
----- Original Message -----
> From: "William Penny" <[log in to unmask]>
> To: "carsten allefeld" <[log in to unmask]>
> Cc: [log in to unmask]
> Sent: Thursday, 6 December, 2018 4:11:06 PM
> Subject: Re: [SPM] Strange AR processes from Bayesian estimation
>
>
>
>
>
> Dear Carsten,
>
>
> Many of these issues are dealt with in the original paper:
>
>
> https://doi.org/10.1016/S1053-8119(03)00071-5
>
>
> For example, section 5.2 of the paper compares (VB versus OLS
> solutions) estimating experimental effect sizes
> by either (I) accounting for autocorrelation structure (the VB
> algorithm) versus (2) ignoring it (OLS). The conclusion being that
> effect sizes are estimated
> more accurately if you also model the autocorrelation of the errors -
> this is the main point of AR modelling.
>
>
> Section 5.3 of the paper tries AR model orders from 0 to 5 and we
> found that an order of 3 was sufficient for all voxels - see map of
> optimal order as a function of location in Figure 8b. This was only
> a single data set however. You can try this on your own data. This
> methodology has been applied in recent papers (fast TR) - see e.g:
>
>
> https://doi.org/10.1016/j.neuroimage.2017.10.043
>
>
> In terms of the algorithm itself, the voxel-wise joint posterior
> distribution over regression coefficients, w, AR coefficients, a,
> and observation noise, lambda, is assumed to have a factorised
> structure: p(w,a,lambda) = q(w)q(a)q(lambda) and the means and
> variances of these approximate posteriors are updated iteratively to
> maximise a lower bound on the model evidence (as with all VB
> approaches).
>
>
> To answer another question, yes, a posterior distribution over the
> regression coefficients is computed, q(a). So the precision of the
> AR estimates is computed. However, these precisions are not stored
> if you run Bayesian 1st level from the GUI (this calls the
> spm_spm_vb.m algorithm which in turn calls spm_vb_glmar (line 838) -
> you can call this latter function directly and analyse small blocks
> of data, looking at optimal model orders, AR precision etc.)
>
>
> All the best,
>
>
> Will.
>
>
>
>
>
>
>
>
>
>
> On Tue, Dec 4, 2018 at 4:12 PM Carsten Allefeld <
> [log in to unmask] > wrote:
>
>
> Dear Will,
>
> thanks for the information!
>
> > In this part of SPM the AR coefficients are estimated using a
> > regression approach in which future time series values (here - the
> > error time series) are regressed onto previous ones. An alternative
> > approach e.g. Burg/Yule-Walker method, first computes the
> > autocorrelation function of the time series and then uses the
> > Yule-Walker relations to estimate the AR coefficients.
>
> Is the regression estimation also implemented in a Bayesian way? If
> yes, is there a way to obtain the posterior precision of the
> estimates, and is any kind of (e.g. spatial smoothness) prior
> implemented? If not, what is the role of AR estimation in "Bayesian
> 1st-level"? Also, it sounds like the regression is applied to
> residuals – how do you make sure the deviation between error
> variance–covariance and residual variance–covariance doesn't bias
> your estimate (what's usually dealt with by ReML)? Sorry, lots of
> questions...
>
> > Perhaps its the case that with higher order models (e.g. 6th order)
> > there are more degrees of freedom to get less well-behaved
> > estimates
> > of the autocorrelation function (sorry to be vague). So it may be
> > worth trying a lower order model (or trying some model order
> > selection here - you could use spm_ar.m from the spectral toolbox
> > to
> > look at other model orders (on the GLM residuals) - with the .fm
> > field returning an approximation to the model evidence).
>
> I only looked at "Bayesian 1st-level" to have voxel-wise estimates of
> autocorrelation, which isn't available in "Classical 1st-level".
> Outside of SPM, the two approaches seem to be AR(6) and ARMA(1,1),
> that's why I selected order 6. Do you recommend order 3?
>
> Best,
> Carsten
>
>
>
> > From: Carsten Allefeld < [log in to unmask] >
> > Sent: 03 December 2018 18:13
> > To: Guillaume Flandin
> > Cc: [log in to unmask] ; William Penny (PSY - Staff)
> > Subject: Re: [SPM] Strange AR processes from Bayesian estimation
> >
> >
> > Dear Guillaume & Will,
> >
> > I noticed a minor inaccuracy in my code transforming AR
> > coefficients
> > into autocorrelation functions. According to comparison with the
> > analytic solution for specific coefficient values, I the corrected
> > version has an absolute error on the order of 1e-16. The updated
> > results have minor changes, but the basic observations are the
> > same.
> >
> > There are 13256 in-mask voxels. Of these, 16 voxels have estimated
> > AR
> > coefficients whose roots' magnitudes are not all below 1
> > (describing
> > nonstationary AR processes).
> >
> > The autocorrelation functions of the remaining 13240 voxels are
> > shown
> > superimposed in the attached autocorrelation.png. As you notice, a
> > non-neglibile number of them have autocorrelation functions that do
> > not decay to 0 over the course of 128 scans (= 256 s).
> >
> > autocorrelation_slice.png shows the autocorrelation over lags 0 to
> > 127 in slice 13 (of 23), containing 1150 in-mask voxels. It is
> > apparent that the voxels with long-range autocorrelation lie mainly
> > at the edge of the brain mask.
> >
> > Guillaume, you asked for time series. timeseries.png shows the
> > underlying BOLD measurements, in the top panel for the 16 voxels
> > with nonstationary AR, and in the lower panel for the 71 voxels
> > where the autocorrelation at lag 127 is larger than 0.1. I can't
> > say
> > that I notice anything special or common to these timeseries.
> >
> > For the moment I will proceed with my data analysis after excluding
> > the 16 + 71 = 87 "weird" voxels. But I think it would be useful if
> > you could look into this further. In particular, in my opinion an
> > AR
> > estimation method should not produce coefficients describing a
> > nonstationary process. Do you agree?
> >
> > Thank you!
> >
> > Best,
> > Carsten
> >
> >
> > ----- Original Message -----
> > > From: "Guillaume Flandin" < [log in to unmask] >
> > > To: "Carsten Allefeld" < [log in to unmask] >
> > > Cc: [log in to unmask] , "William Penny (PSY)" <
> > > [log in to unmask] >
> > > Sent: Wednesday, 28 November, 2018 11:35:53 AM
> > > Subject: Re: [SPM] Strange AR processes from Bayesian estimation
> > >
> > > Dear Carsten,
> > >
> > > All the voxels exhibiting pathological behaviour seem to be at
> > > the
> > > boundary of the brain mask so, as you say, simply discarding them
> > > is
> > > probably the easiest thing to do. To understand what is going on,
> > > it
> > > would be useful to extract and display their time series.
> > > I copy this email to Will so that he can add any further comments
> > > or
> > > advice on the spatial noise prior.
> > >
> > > Best regards,
> > > Guillaume.
> > >
> > >
> > > On 27/11/2018 19:20, Carsten Allefeld wrote:
> > > > Dear Guillaume,
> > > >
> > > > thanks for replying!
> > > >
> > > >> These results were obtained with which spatial noise prior
> > > >> option?
> > > >> The
> > > >> interface lists five options: UGL, GMRF, LORETA, Tissue-type
> > > >> and
> > > >> Robust.
> > > >
> > > > I used the default, "UGL". Which one would you recommend?
> > > >
> > > > Options in detail:
> > > > fmri_est.method.Bayesian.space.volume.block_type = 'Slices';
> > > > fmri_est.method.Bayesian.signal = 'UGL';
> > > > fmri_est.method.Bayesian.ARP = 6;
> > > > fmri_est.method.Bayesian.noise.UGL = 1;
> > > > fmri_est.method.Bayesian.LogEv = 'No';
> > > > fmri_est.method.Bayesian.anova.first = 'No';
> > > > fmri_est.method.Bayesian.anova.second = 'No';
> > > > fmri_est.method.Bayesian.gcon = struct('name', {}, 'convec',
> > > > {});
> > > >
> > > >> Could you show a map of where the voxels you are concerned
> > > >> about
> > > >> are?
> > > >
> > > > Attached are plots of the absolute value of the autocorrelation
> > > > at
> > > > lags 0 to 127 in a middle slice (#13).
> > > > Comparison with the coregistered T1 indicates that they are
> > > > located
> > > > mainly at the outer edge of gray matter (maybe meninges), but
> > > > also
> > > > frontally and posteriorly slightly into the longitudinal
> > > > fissure.
> > > >
> > > > That suggests I should simply exclude these voxels from further
> > > > analysis.
> > > >
> > > > Do you have a suggestion which criterion to use?
> > > > The data-based threshold (2.02e-9) discards more than half of
> > > > the
> > > > brain, and any other threshold seems arbitrary.
> > > >
> > > > Best,
> > > > Carsten
> > > >
> > > >
> > > >>
> > > >> Best regards,
> > > >> Guillaume.
> > > >>
> > > >>
> > > >> On 27/11/2018 17:11, Carsten Allefeld wrote:
> > > >>> Hello all,
> > > >>>
> > > >>> I'm interested in getting local estimates of temporal
> > > >>> autocorrelation in SPM, and for that purpose used Bayesian
> > > >>> 1st-level estimation.
> > > >>> The fMRI data I used to test that have 3 sessions of 128
> > > >>> scans
> > > >>> at
> > > >>> a
> > > >>> TR of 2 s and 64x64x23 voxels of size 4x4x4 mm, unsmoothed,
> > > >>> of
> > > >>> which approximately 13,000 are within brain.
> > > >>>
> > > >>> I then extracted the AR coefficients (order 6) for the first
> > > >>> session (Sess1_AR_0001.nii to Sess1_AR_0006.nii) and used the
> > > >>> Yule–Walker equations iteratively to obtain the corresponding
> > > >>> autocorrelation function across lags 0 to 128.
> > > >>>
> > > >>> The results are strange (see attached plot):
> > > >>> – In 16 voxels the AR coefficients describe a non-stationary
> > > >>> process. After excluding them:
> > > >>> – At lag 127, 7906 voxels have an autocorrelation > 1e-6,
> > > >>> 1651
> > > >>> voxels > 1e-3, and 113 voxels > 0.1.
> > > >>> – The largest negative autocorrelation at lag 127 is
> > > >>> -2.02e-9.
> > > >>> If
> > > >>> I
> > > >>> take that as an indicator of numerical/estimation precision,
> > > >>> there
> > > >>> are 8185 voxels where the autocorrelation at lag 127 is
> > > >>> different
> > > >>> from 0 (> +2.02e-9).
> > > >>>
> > > >>> This makes me suspect that the AR estimation in "Bayesian
> > > >>> 1st-level" is not very reliable. Is there something I might
> > > >>> have
> > > >>> done wrong?
> > > >>>
> > > >>> Is there a recommended postprocessing for the AR coefficients
> > > >>> or
> > > >>> autocorrelation functions?
> > > >>> I thought about tapering à la FSL, or clustering as a crude
> > > >>> form
> > > >>> of
> > > >>> spatial regularization.
> > > >>> Or should I simply exclude voxels with unbelievably
> > > >>> long-range
> > > >>> autocorrelation?
> > > >>>
> > > >>> Thank!
> > > >>>
> > > >>> Best,
> > > >>> Carsten
> > > >>>
> > > >>
> > > >> --
> > > >> Guillaume Flandin, PhD
> > > >> Wellcome Centre for Human Neuroimaging
> > > >> UCL Queen Square Institute of Neurology
> > > >> London WC1N 3BG
> > > >>
> > >
> > > --
> > > Guillaume Flandin, PhD
> > > Wellcome Centre for Human Neuroimaging
> > > UCL Queen Square Institute of Neurology
> > > London WC1N 3BG
> > >
> >
>
|