FAQ - OHBA-analysis/HMM-MAR GitHub Wiki

Slowly growing FAQ – will keep adding question as they come up.

How do I get started?

We have applied the HMM to various types of data, including M/EEG and fMRI. For M/EEG, a good point to start is Quinn et al. (2018), and the associated tutorial. More examples can be found in the Examples section of this Wiki, and here – which can be easily tailored to other data sets and settings. For fMRI, this one can be used as a starting point.

Can I get in contact and ask questions?

Yes, please. But maybe after having a look to this FAQ and the rest of the Wiki.

How do I choose the type of model?

That depends on the type of data and what aspects of the data one wants to focus on. If analysing MEG or EEG data with many sensors or sources, then the Time-Delay Embedded HMM is a sensible choice (options.embeddedlags) If using intracranial data or just a few sensors or channels, and especially if you care about directionality, then the HMM-MAR is adequate; note that this is a more complex model, so, if that causes problems, the Time-Delay Embedded HMM is a safe bet. When analysing fMRI data, we have always used an HMM with Gaussian-distributed states. If one just cares about connectivity, then the mean of the Gaussian distribution can be fixed to zero (options.zeromean=1). If the registration is not very good, then maybe the covariance matrix can be set to be equal for all states (options.covtype=’uniquefull’); see also Q12.

If the code is able to use different state distributions, then why is it called HMM-MAR?

The MAR (Multivariate AutoRegressive model) was the model used in our first publication with this toolbox, so in part this is for historic reasons. Also, the rest of the distributions can be considered a generalisation of the MAR.

Why, every time I run it, I don’t get exactly the same solution?

Same than with ICA and many other computational methods, the inference process that is used to estimate the parameters of the model can be a bit complex. The optimisation problem that this process is trying to solve often has different local minima, and, since there is some randomness in the initial conditions of the inference, this can converge differently in each run. This depends on the data and the specification of the model, but, in our experience, the optimisation problem often has one single solution, sometimes two, and only occasionally more than that. In any case, these solutions are always valid, the information is there, but is represented slightly differently. For this reason, we would recommend run the model various times and inspect every solution separately.

We can also use this stochasticity to our advantage. For example, if we want to find differences in certain feature (e.g. fractional occupancy) between two groups, we could combine the different runs into a single test to boost statistical power; see Vidaurre et al., 2019 for details about how to to that in practice, and this Github repository

It takes too long to run, what can I do?

Using stochastic inference will make things much faster; see Vidaurre et al. 2018 and check the corresponding section of the User Guide.

And if it uses too much memory?

Same answer.

How many states should I choose?

There is not really an exact answer to this one. More states will give a more fine-grained solution; fewer states will be coarser. It depends on the question and data, but, in general, running the model in half-splits of the data and comparing the robustness of the estimations is a good way to decide an upper bound in the number of states.

Which parcellation should I use in MEG? How many parcels? Should I use subcortical ones?

When looking at whole brain, it is normal to use a 30-60 regions parcellation. This recommendation is not specific to the HMM, but more due to the fact that we normally use leakage correction (correction for volume conduction), which removes a lot of variance from the signal. A higher number will result in removing too much signal for an HMM analysis (or any other analysis) to be useful. Of course, one should always use a parcellation that suits the scientific question at hand.

About the subcortical areas, the signal-to-noise ratio for those is obviously lower, so it would make sense to check the SNR first. If that is too low, then adding the subcortical areas would just make the HMM estimation noisier – but still primarily driven by the regions with a higher SNR (e.g. motor)

How can I interrogate the properties of the states?

Depending on the type of model used, there are various functions that can be readily invoked. If the model is a Gaussian distribution endowed with a mean parameter (the model was trained with options.zeromean=0), then the function getMean will extract the mean activation for any given state.

If the model is a Gaussian distribution endowed a state-specific covariance matrix (the model was trained with options.covtype=full, or options.covtype=diag), then the function getFuncConn will do the job.

If the model is a MAR model (i.e. options.order was higher than 0), then the function getMARmodel can be used.

If the model is a time-delay embedded HMM (TDE-HMM; the model was trained specifying options.embeddedlags), then getAutoCovMat will provide the autocovariance (or cross-correlation, depending on terminology) matrix of the state in question.

Do you really think the brain has a discrete number of states?

I don’t know. This is just an analysis technique, giving you a certain perspective of the data. Whether the HMM is useful or not, depends on the aspect in the data that we aim to interrogate.

I am getting precision problems, what can I do?

Sometimes warnings saying things like “Matrix is close to singular or badly scaled”, or an error informing of “Issues of precision”, can occur. There are two reasons because this might happen.

First, it is possible that your data has some huge artefact, like for example a large peak. If that is the case, then the likelihood for the states to explain this part of the data is so small that falls under the threshold of what Matlab can represent. When this happens, one normally sees the “singular matrix” warning. The only solution in this case is to clean these artefacts.

Second, it is possible that you are using a too complex model (e.g. a MAR model with too large order), which explains the data so well that the error falls under the precision threshold. This happens also when we are using a MAR model on too many channels. The solution is either to reduce the order of the MAR, or to use PCA (options.pca) on the channels, or to use the Time-Delay Embedded HMM.

One state is dominating the entire data set.

The chosen model is too complex, and one single state explain a very large amount of variance so that the others become unnecessary. Most likely you are using the MAR model. Again: reduce the order of the MAR, or to use PCA (options.pca) on the channels, or to use the Time-Delay Embedded HMM.

Single states explain entire subjects or sessions, so that there is little or none state switching.

This typically happens because the registration of the subjects is not very good, or/and the amount of data per subject is not large enough. Specially, when using full covariance matrices (covtype=’full’), it might be that the (regions by regions) average functional connectivity matrices are so dissimilar between subjects that the best solution the HMM can come up with is to assign one state to each subject; that is, between-subject differences massively dominate within-subject changes. The function corrStaticFC can help to diagnose if this is the case (see User Guide).

In fMRI, the best solution is to improve the registration, or, if not possible, to use covtype=’uniquefull’, which will make all the states to have the same connectivity; see User Guide. Another option to keep the functional connectivity information is to use PCA (options.pca), which we have found to help sometimes. We are actually working to provide a better solution for this problem, but that's still work in progress.

When using models that contain spectral information (MAR or Time-Delay Embedded) this might also happen when the frequency profiles of the subjects are too distinct – but this is rare. Another possibility, when working with source-reconstructed M/EEG data, is that the sign ambiguity issue is causing the problem (see the "Sign flipping" section in the User Guide).

I am getting errors on the parfor loops (parallel computing), what is going on?

Matlab's parallel computing can be a temperamental creature, and also the error messages are not great to diagnose what happens. We have found errors for certain combinations of Matlab distributions and operating systems, which do not seem to be related to the HMM code. If the problem persists, disabling the parallel computing feature by setting options.useParallel=0 should solve the issue.

In fMRI, how can I focus my decomposition on functional connectivity?

Changes in the variance of the signal are always going to affect the estimation, but the best we can do (for now) is to use options.zeromean=1 and options.covtype=’full’.

Can I control the switching rate of the states?

Yes, with options.DirichletDiag, which tweaks the prior probability of remaining on the same state. Unfortunately, this depends on how much data we have (how many time points), so it is difficult to provide a solid recipe. You might want to increase one order of magnitude at a time, starting from options.DirichletDiag=10.

The length of the state time courses does not coincide with the data, what can I do?

This happens when using the MAR model or the Time-Delay Embedded HMM. The function padGamma can fill the state time courses to have the same length of the data, for convenience.