spm_mvb_ui.m 8.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235
  1. function [MVB] = spm_mvb_ui(xSPM,SPM,MVB)
  2. % Multivariate Bayes (Bayesian decoding of a contrast)
  3. % FORMAT [MVB] = spm_mvb_ui(xSPM,SPM,MVB)
  4. %
  5. % Sets up, evaluates and saves an MVB structure:
  6. %
  7. % MVB.contrast % contrast structure
  8. % MVB.name % name
  9. % MVB.c % contrast weight vector
  10. % MVB.M % MVB model (see below)
  11. % MVB.X % subspace of design matrix
  12. % MVB.Y % multivariate response
  13. % MVB.X0 % null space of design
  14. % MVB.XYZ % location of voxels (mm)
  15. % MVB.V % serial correlation in response
  16. % MVB.K % whitening matrix
  17. % MVB.VOX % voxel scaling
  18. % MVB.xyzmm % centre of VOI (mm)
  19. % MVB.Space % VOI definition
  20. % MVB.Sp_info % parameters of VOI
  21. % MVB.Ni % number of greedy search steps
  22. % MVB.sg % size of reedy search split
  23. % MVB.priors % model (spatial prior)
  24. % MVB.fSPM % SPM analysis (.mat file)
  25. %
  26. % where MVB.M contains the following fields:
  27. %
  28. % F: log-evidence [F(0), F(1),...]
  29. % G: covariance partition indices
  30. % h: covariance hyperparameters
  31. % U: ordered patterns
  32. % qE: conditional expectation of voxel weights
  33. % qC: conditional variance of voxel weights
  34. % Cp: prior covariance (ordered pattern space)
  35. % cp: prior covariance (original pattern space)
  36. %
  37. %--------------------------------------------------------------------------
  38. % This routine uses a multivariate Bayesian (MVB) scheme to decode or
  39. % recognise brain states from neuroimages. It resolves the ill-posed
  40. % many-to-one mapping, from voxel values or data features to a target
  41. % variable, using a parametric empirical or hierarchical Bayesian model.
  42. % This model is inverted using standard variational techniques, in this
  43. % case expectation maximisation, to furnish the model evidence and the
  44. % conditional density of the model's parameters. This allows one to compare
  45. % different models or hypotheses about the mapping from functional or
  46. % structural anatomy to perceptual and behavioural consequences (or their
  47. % deficits). The aim of MVB is not to predict (because the outcomes are
  48. % known) but to enable inference on different models of structure-function
  49. % mappings; such as distributed and sparse representations. This allows one
  50. % to optimise the model itself and produce predictions that outperform
  51. % standard pattern classification approaches, like support vector machines.
  52. % Technically, the model inversion and inference uses the same empirical
  53. % Bayesian procedures developed for ill-posed inverse problems (e.g.,
  54. % source reconstruction in EEG).
  55. %
  56. % CAUTION: MVB should not be used to establish a significant mapping
  57. % between brain states and some classification or contrast vector. Its use
  58. % is limited to comparison of different models under the assumption
  59. % (hyperprior) that this mapping exists. To ensure the mapping exists, use
  60. % CVA or compute the randomisation p-value (see spm_mvb_p)
  61. %
  62. % See: spm_mvb and
  63. %
  64. % Bayesian decoding of brain images.
  65. % Friston K, Chu C, Mourao-Miranda J, Hulme O, Rees G, Penny W, Ashburner J.
  66. % Neuroimage. 2008 Jan 1;39(1):181-205
  67. %
  68. % Multiple sparse priors for the M/EEG inverse problem.
  69. % Friston K, Harrison L, Daunizeau J, Kiebel S, Phillips C, Trujillo-Barreto
  70. % N, Henson R, Flandin G, Mattout J.
  71. % Neuroimage. 2008 Feb 1;39(3):1104-20.
  72. %
  73. % Characterizing dynamic brain responses with fMRI: a multivariate approach.
  74. % Friston KJ, Frith CD, Frackowiak RS, Turner R.
  75. % Neuroimage. 1995 Jun;2(2):166-72.
  76. %__________________________________________________________________________
  77. % Copyright (C) 2007-2017 Wellcome Trust Centre for Neuroimaging
  78. % Karl Friston
  79. % $Id: spm_mvb_ui.m 7162 2017-08-30 11:47:07Z guillaume $
  80. %-Get figure handles and set title
  81. %--------------------------------------------------------------------------
  82. Finter = spm_figure('FindWin','Interactive');
  83. spm_results_ui('Clear');
  84. spm_input('!DeleteInputObj');
  85. header = get(Finter,'Name');
  86. spm_clf(spm_figure('FindWin','MVB'));
  87. %-Get contrast: only the first line of F-contrast
  88. %--------------------------------------------------------------------------
  89. try
  90. contrast = SPM.xCon(xSPM.Ic).name;
  91. c = SPM.xCon(xSPM.Ic).c(:,1);
  92. catch
  93. contrast = MVB.contrast;
  94. c = MVB.c;
  95. end
  96. %-Get VOI name
  97. %--------------------------------------------------------------------------
  98. try
  99. name = MVB.name;
  100. catch
  101. name = spm_input('name','-8','s',contrast);
  102. end
  103. name = strrep(name,' ','_');
  104. name = ['MVB_' name];
  105. %-Get current location {mm}
  106. %--------------------------------------------------------------------------
  107. try
  108. xyzmm = MVB.xyzmm;
  109. catch
  110. xyzmm = spm_results_ui('GetCoords');
  111. end
  112. %-Specify search volume
  113. %--------------------------------------------------------------------------
  114. try
  115. xY = MVB.xY;
  116. MVB = rmfield(MVB,'xY');
  117. catch
  118. xY = [];
  119. end
  120. xY.xyz = xyzmm;
  121. Q = ones(1,size(SPM.xVol.XYZ,2));
  122. XYZmm = SPM.xVol.M(1:3,:)*[SPM.xVol.XYZ; Q];
  123. [xY,XYZ,j] = spm_ROI(xY, XYZmm);
  124. %-Get explanatory variables (data)
  125. %--------------------------------------------------------------------------
  126. XYZ = XYZmm(:,j);
  127. Y = spm_get_data(SPM.xY.VY,SPM.xVol.XYZ(:,j));
  128. %-Check there are intracranial voxels
  129. %--------------------------------------------------------------------------
  130. if isempty(Y)
  131. spm('alert*',{'No voxels in this VOI';'Please use a larger volume'},...
  132. 'Multivariate Bayes');
  133. end
  134. %-Get model[s]
  135. %--------------------------------------------------------------------------
  136. try
  137. priors = lower(MVB.priors);
  138. catch
  139. str = {'compact','sparse','smooth','support'};
  140. Ip = spm_input('model (spatial prior)','!+1','m',str);
  141. priors = str{Ip};
  142. end
  143. %-Number of iterations
  144. %--------------------------------------------------------------------------
  145. try
  146. sg = MVB.sg;
  147. catch
  148. str = 'size of successive subdivisions';
  149. sg = spm_input(str,'!+1','e',.5);
  150. end
  151. %-MVB is now specified
  152. %==========================================================================
  153. spm('Pointer','Watch')
  154. %-Get target and confounds
  155. %--------------------------------------------------------------------------
  156. X = SPM.xX.X;
  157. X0 = X*(speye(length(c)) - c*pinv(c));
  158. try
  159. % accounting for multiple sessions
  160. %----------------------------------------------------------------------
  161. tmpX0 = [];
  162. for ii = 1:length(SPM.xX.K)
  163. tmp = zeros(sum(SPM.nscan),size(SPM.xX.K(ii).X0,2));
  164. tmp(SPM.xX.K(ii).row,:) = SPM.xX.K(ii).X0;
  165. tmpX0 = [tmpX0 tmp];
  166. end
  167. X0 = [X0 tmpX0];
  168. end
  169. X = X*c;
  170. %-Serial correlations
  171. %--------------------------------------------------------------------------
  172. V = SPM.xVi.V;
  173. %-Invert
  174. %==========================================================================
  175. VOX = diag(abs(SPM.xVol.M));
  176. U = spm_mvb_U(Y,priors,X0,XYZ,VOX);
  177. try
  178. Ni = MVB.Ni;
  179. catch
  180. str = 'Greedy search steps';
  181. Ni = spm_input(str,'!+1','i',max(8,ceil(log(size(U,2))/log(1/sg))));
  182. end
  183. M = spm_mvb(X,Y,X0,U,V,Ni,sg);
  184. M.priors = priors;
  185. %-Assemble results
  186. %--------------------------------------------------------------------------
  187. MVB.contrast = contrast; % contrast of interest
  188. MVB.name = name; % name
  189. MVB.c = c; % contrast weight vector
  190. MVB.M = M; % MVB model (see below)
  191. MVB.X = X; % subspace of design matrix
  192. MVB.Y = Y; % multivariate response
  193. MVB.X0 = X0; % null space of design
  194. MVB.XYZ = XYZ; % location of voxels (mm)
  195. MVB.V = V; % serial correlation in repeosne
  196. MVB.K = full(V)^(-1/2); % whitening matrix
  197. MVB.VOX = SPM.xVol.M; % voxel scaling
  198. MVB.xyzmm = xyzmm; % centre of VOI (mm)
  199. MVB.Space = xY.def; % VOI definition
  200. MVB.Sp_info = xY.spec; % parameters of VOI
  201. MVB.Ni = Ni; % number of greedy search steps
  202. MVB.sg = sg; % size of reedy search split
  203. MVB.priors = priors; % model (spatial prior)
  204. MVB.fSPM = fullfile(SPM.swd,'SPM.mat'); % SPM analysis (.mat file)
  205. %-Display
  206. %==========================================================================
  207. if ~spm('CmdLine'), spm_mvb_display(MVB); end
  208. %-Save
  209. %--------------------------------------------------------------------------
  210. save(fullfile(SPM.swd,[name '.mat']),'MVB', spm_get_defaults('mat.format'))
  211. assignin('base','MVB',MVB)
  212. %-Reset title
  213. %--------------------------------------------------------------------------
  214. set(Finter,'Name',header)
  215. spm('Pointer','Arrow')