Homec4science
Diffusion Slepians (master) Browse

Slepi/Code/Utilities/spm12/spm_eeg_invert.m
8a186cd82567master

Last Change
Raw File

spm_eeg_invert.m
View Options

function [D] = spm_eeg_invert(D, val)
% ReML inversion of multiple forward models for EEG-MEG
% FORMAT [D] = spm_eeg_invert(D)
% ReML estimation of regularisation hyperparameters using the
% spatiotemporal hierarchy implicit in EEG/MEG data
%
% Requires:
% D{i}.inv{val}.inverse:
%
% inverse.modality - modality to use in case of multimodal datasets
%
% inverse.trials - D.events.types to invert
% inverse.type - 'GS' Greedy search on MSPs
% 'ARD' ARD search on MSPs
% 'MSP' GS and ARD multiple sparse priors
% 'LOR' LORETA-like model
% 'IID' minimum norm
% inverse.woi - time window of interest ([start stop] in ms)
% inverse.lpf - band-pass filter - low frequency cut-off (Hz)
% inverse.hpf - band-pass filter - high frequency cut-off (Hz)
% inverse.Han - switch for Hanning window
% inverse.xyz - (n x 3) locations of spherical VOIs
% inverse.rad - radius (mm) of VOIs
%
% inverse.Nm - maximum number of channel modes
% inverse.Nr - maximum number of temporal modes
% inverse.Np - number of sparse priors per hemisphere
% inverse.smooth - smoothness of source priors (0 to 1)
% inverse.Na - number of most energetic dipoles
% inverse.sdv - standard deviations of Gaussian temporal correlation
% inverse.pQ - any source priors (e.g. from fMRI); vector or matrix
% inverse.Qe - any sensor error components (e.g. empty-room data)
% inverse.dplot - make diagnostics plots (0 or 1)
% inverse.STAT - flag for stationarity assumption, which invokes a
% full DCT temporal projector (from lpf to hpf Hz)
%
% Evaluates:
%
% inverse.M - MAP projector (reduced)
% inverse.J{i} - Conditional expectation (i conditions) J = M*U*Y
% inverse.L - Lead field (reduced UL := U*L)
% inverse.qC - spatial covariance
% inverse.qV - temporal correlations
% inverse.T - temporal projector
% inverse.U(j) - spatial projector (j modalities)
% inverse.Y{i} - reduced data (i conditions) UY = UL*J + UE
% inverse.Is - Indices of active dipoles
% inverse.It - Indices of time bins
% inverse.Ic{j} - Indices of good channels (j modalities)
% inverse.Nd - number of dipoles
% inverse.pst - peristimulus time
% inverse.dct - frequency range
% inverse.F - log-evidence
% inverse.VE - variance explained in spatial/temporal subspaces (%)
% inverse.R2 - variance in subspaces accounted for by model (%)
% inverse.scale - scaling of data for each of j modalities
%__________________________________________________________________________
%
% 1. This routine implements "group-based" inversion, corresponding to
% ill-posed linear models of the following form:
%
% [AY{1}...AY{n}] = L(1} * [J{1}...J{n}] + [e{1}...e{n}]
%
% where AY{i} are the spatially normalized or adjusted data from subject i
% that would have been seen if the lead-field L{i} = L{1}. The ensuing
% Gaussian process priors on sources are then used to estimate subject-
% specific MAP estimates of J{i} using
%
% AY{i} = L(1} * J{i} + e{i}
%
% using spatial priors from the group model above.
%
% Here, A{i} = L{1}*pinv(L{i}) =>
% AY{i} = A(i}*L(i}*J{i}
% = L(1}*J{i}
%
% Potential scaling differences between the lead-fields are handled by
% scaling L{1} such that trace(L{1}*L{1}') = constant (number of spatial
% modes or channels), while scaling the data such that trace(AY{n}*AY{n}') =
% constant over subjects (and modalities; see below).
%
% See: Electromagnetic source reconstruction for group studies.
% Litvak V, Friston K.
% NeuroImage. 2008 Oct 1;42(4):1490-8.
%
%__________________________________________________________________________
%
% 2. It also implements "fusion" of different types of MEG and EEG data,
% corresponding to ill-posed linear models of the following form:
%
% AY{1}{1,...,t} = L(1} * J{1,...,t} + e{{1,...,t}}
% AY{2}{1,...,t} = L(2} e{{2,...,t}}
% .
% .
% .
% AY{m}{1,...,t} = L(n} e{{n,...,t}}
%
% Under empirical priors on J{1,...,t} for m modalities with t trial types.
%
% See: MEG and EEG data fusion: Simultaneous localisation of face-evoked
% responses.
% Henson R, Mouchlianitis E & Friston K.
% Neuroimage. 2009. 47:581-9.
%__________________________________________________________________________
%
% 3. It also allows incorporation of spatial source priors, eg, from fMRI
% (see spm_eeg_inv_fmripriors.m). Note that if a vector is passed in
% inverse.pQ, then variance components used (pass a matrix if a covariance
% component is desired).
%
% See: A Parametric Empirical Bayesian framework for fMRI-constrained
% MEG/EEG source reconstruction.
% Henson R, Flandin G, Friston K & Mattout J.
% Human Brain Mapping (in press).
%__________________________________________________________________________
%
% The routine essentially consists of two steps:
%
% 1. Optimisation of spatial source priors over subjects
% 2. Re-inversion of each subject, fusing across all modalities
%__________________________________________________________________________
% Copyright (C) 2006-2014 Wellcome Trust Centre for Neuroimaging
% Karl Friston
% $Id: spm_eeg_invert.m 6182 2014-09-18 12:03:18Z guillaume $
% check whether this is a group inversion for (Nl) number of subjects
%--------------------------------------------------------------------------
if ~iscell(D), D = {D}; end
Nl = length(D);
% D - SPM data structure
%==========================================================================
if nargin > 1
D{1}.val = val;
elseif ~isfield(D{1}, 'val')
D{1}.val = 1;
end
for i = 2:Nl
D{i}.val = D{1}.val;
end
% forward model
%--------------------------------------------------------------------------
inverse = D{1}.inv{D{1}.val}.inverse;
% defaults
%--------------------------------------------------------------------------
try, STAT = inverse.STAT; catch, STAT = 0; end
try, type = inverse.type; catch, type = 'GS'; end
try, s = inverse.smooth; catch, s = 0.6; end
try, Np = inverse.Np; catch, Np = 256; end
try, Nr = inverse.Nr; catch, Nr = 16; end
try, xyz = inverse.xyz; catch, xyz = [0 0 0]; end
try, rad = inverse.rad; catch, rad = 128; end
try, lpf = inverse.lpf; catch, lpf = 0; end
try, hpf = inverse.hpf; catch, hpf = 48; end
try, sdv = inverse.sdv; catch, sdv = 4; end
try, Han = inverse.Han; catch, Han = 1; end
try, woi = inverse.woi; catch, woi = []; end
try, pQ = inverse.pQ; catch, pQ = []; end
try, dp = inverse.dplot; catch, dp = 0; end
% get specified modalities to invert (default to all)
%--------------------------------------------------------------------------
try
modalities = inverse.modality;
if ~iscell(modalities)
modalities = {modalities};
end
catch
for m = 1:length(D{1}.inv{D{1}.val}.forward)
modalities{m} = D{1}.inv{D{1}.val}.forward(m).modality;
end
end
Nmod = numel(modalities); % number of modalities
Nmax = Nr; % max number of temporal modes
% check lead fields and get number of dipoles (Nd) and channels (Nc)
%==========================================================================
for i = 1:Nl
fprintf('Checking lead fields for subject %i\n',i)
[L,D{i}] = spm_eeg_lgainmat(D{i});
for m = 1:Nmod
% Check gain or lead-field matrices
%------------------------------------------------------------------
Ic{i,m} = indchantype(D{i}, modalities{m}, 'GOOD');
Nd(i) = size(L,2);
Nc(i,m) = length(Ic{i,m});
if isempty(Ic{i,m})
errordlg(['Modality ' modalities{m} 'is missing from file ' D{i}.fname]);
return
end
if any(diff(Nd))
errordlg('Please ensure subjects have the same number of dipoles')
return
end
% Check for null space over sensors (SX) and remove it
%------------------------------------------------------------------
try
SX = D{i}.sconfounds{m};
R{i,m} = speye(Nc(i,m),Nc(i,m)) - SX*spm_pinv(SX);
catch
R{i,m} = speye(Nc(i,m),Nc(i,m));
end
end
end
fprintf(' - done\n')
% Compute spatial coherence: Diffusion on a normalised graph Laplacian GL
%==========================================================================
fprintf('Computing Green function from graph Laplacian:')
%--------------------------------------------------------------------------
Nd = Nd(1); % number of dipoles
vert = D{1}.inv{D{1}.val}.mesh.tess_mni.vert;
face = D{1}.inv{D{1}.val}.mesh.tess_mni.face;
A = spm_mesh_distmtx(struct('vertices',vert,'faces',face),0);
GL = A - spdiags(sum(A,2),0,Nd,Nd);
GL = GL*s/2;
Qi = speye(Nd,Nd);
QG = sparse(Nd,Nd);
for i = 1:8
QG = QG + Qi;
Qi = Qi*GL/i;
end
QG = QG.*(QG > exp(-8));
QG = QG*QG;
clear Qi A GL
fprintf(' - done\n')
% check for (e.g., empty-room) sensor components (in Qe{1})
%==========================================================================
QE = cell(Nl,Nmod);
for i = 1:Nl
for m = 1:Nmod
try
QE{i,m} = D{i}.inv{D{i}.val}.inverse.Qe{m};
QE{i,m} = Nc(i,m)*QE{i,m}/trace(QE{i,m});
if length(QE{i,m}) ~= Nc(i,m)
errordlg('error component (modality %s; subject %d) does not match number of channels (%d)\n',modalities{m},i,Nc(i,m))
return
end
fprintf('Using sensor error component provided...\n');
% assume i.i.d. if not specified
%------------------------------------------------------------------
catch
QE{i,m} = 1;
end
end
end
%==========================================================================
% Spatial projectors (adjusting for different Lead-fields)
%==========================================================================
fprintf('Optimising and aligning spatial modes ...\n')
% Use recursive (regularised) least squares to find average lead-field
%--------------------------------------------------------------------------
Is = 1:Nd;
UL = cell(Nmod,1);
for m = 1:Nmod
% Initialise average lead-field with L{1}
%----------------------------------------------------------------------
UL{m} = R{1,m}*spm_eeg_lgainmat(D{1},Is,D{1}.chanlabels(Ic{1,m}));
AA = 1;
% pre-compute regularised inverses (for speed)
%----------------------------------------------------------------------
for i = 1:Nl
L = R{i,m}*spm_eeg_lgainmat(D{i},Is,D{i}.chanlabels(Ic{i,m}));
iL{i} = spm_inv(L*L');
end
% Optimise alignment matrices A such that A{m}*L{m} = <A{m}*L{m}>m
%----------------------------------------------------------------------
for j = 1:8
% eliminate redundant virtual channels
%------------------------------------------------------------------
fprintf('Aligning - iteration: %i\n',j)
UL{m} = spm_sqrtm(spm_inv(AA))*UL{m};
% eliminate low SNR spatial modes
%------------------------------------------------------------------
U = spm_svd((UL{m}*UL{m}'));
UL{m} = U'*UL{m};
Nm(m) = size(UL{m},1);
% normalise lead-field
%------------------------------------------------------------------
Scale = sqrt(trace(UL{m}*UL{m}')/Nm(m));
UL{m} = UL{m}/Scale;
% spatial projectors A{i,m) for i = 1,...,Nl subjects
%------------------------------------------------------------------
AL = 0;
AA = 0;
for i = 1:Nl
L = R{i,m}*spm_eeg_lgainmat(D{i},Is,D{i}.chanlabels(Ic{i,m}));
A{i,m} = UL{m}*L'*iL{i};
AL = AL + A{i,m}*L;
AA = AA + A{i,m}*A{i,m}';
end
% re-compute average
%------------------------------------------------------------------
UL{m} = AL/Nl;
if Nl == 1, break, end
end
% Report
%----------------------------------------------------------------------
fprintf('Using %d spatial modes for modality %s\n',Nm(m),modalities{m})
if dp
spm_figure('GetWin','Lead fields');
for i = 1:Nl
L = R{i,m}*spm_eeg_lgainmat(D{i},Is(1:8:end),D{i}.chanlabels(Ic{i,m}));
subplot(Nl,4,(i - 1)*4 + 1)
imagesc(A{i,m}*A{i,m}'), axis square
subplot(Nl,4,(i - 1)*4 + 2)
imagesc(A{i,m}'*A{i,m}), axis square
subplot(Nl,4,(i - 1)*4 + 3)
plot(sum(A{i,m}.^2,2)), axis square
subplot(Nl,4,(i - 1)*4 + 4)
plot(L'*A{i,m}'), axis square tight
end
drawnow
end
end
% check restriction: assume radii are the same for all (Nv) VOI
%==========================================================================
Nv = size(xyz,1);
if length(rad) ~= Nv
rad = rad(1)*ones(Nv,1);
else
rad = rad(:);
end
% Restrict source space to Ns sources by eliminating dipoles
%--------------------------------------------------------------------------
Is = sparse(Nd,1);
for i = 1:Nv
Iv = sum([vert(:,1) - xyz(i,1), ...
vert(:,2) - xyz(i,2), ...
vert(:,3) - xyz(i,3)].^2,2) < rad(i)^2;
Is = Is | Iv;
end
Is = find(Is);
vert = vert(Is,:);
QG = QG(Is,Is);
for m = 1:Nmod
UL{m} = UL{m}(:,Is);
end
Ns = length(Is);
%==========================================================================
% Temporal projector
%==========================================================================
% loop over Nl lead-fields (subjects)
%--------------------------------------------------------------------------
Nn = zeros(1,Nl); % number of samples
AY = {}; % pooled response for MVB
AYYA = 0; % pooled response for ReML
for i = 1:Nl
% Time-window of interest
%----------------------------------------------------------------------
if isempty(woi)
w{i} = 1000*[min(D{i}.time) max(D{i}.time)];
else
w{i} = woi;
end
It{i} = (w{i}/1000 - D{i}.timeonset)*D{i}.fsample + 1;
It{i} = max(1,It{i}(1)):min(It{i}(end), length(D{i}.time));
It{i} = fix(It{i});
% Peristimulus time
%----------------------------------------------------------------------
pst{i} = 1000*D{i}.time; % peristimulus time (ms)
pst{i} = pst{i}(It{i}); % windowed time (ms)
dur = (pst{i}(end) - pst{i}(1))/1000; % duration (s)
dct{i} = (It{i} - It{i}(1))/2/dur; % DCT frequencies (Hz)
Nb(i) = length(It{i}); % number of time bins
% Serial correlations
%----------------------------------------------------------------------
K = exp(-(pst{i} - pst{i}(1)).^2/(2*sdv^2));
K = toeplitz(K);
qV{i} = sparse(K*K');
% Confounds and temporal subspace
%----------------------------------------------------------------------
T = spm_dctmtx(Nb(i),Nb(i));
j = find( (dct{i} >= lpf) & (dct{i} <= hpf) );
T = T(:,j);
dct{i} = dct{i}(j);
% Hanning operator (if requested)
%----------------------------------------------------------------------
if Han
W = sparse(1:Nb(i),1:Nb(i),spm_hanning(Nb(i)));
else
W = 1;
end
% get trials or conditions
%----------------------------------------------------------------------
try
trial = D{i}.inv{D{i}.val}.inverse.trials;
catch
trial = D{i}.condlist;
end
Nt(i) = length(trial);
% get temporal covariance (Y'*Y) to find temporal modes
%======================================================================
MY = cell(Nmod,1); % mean response
YTY = sparse(0); % accumulator
for m = 1:Nmod % loop over modalities
% get (spatially aligned) data
%------------------------------------------------------------------
N = 0;
YY = 0;
MY{m} = 0;
for j = 1:Nt(i) % pool over conditions
c = D{i}.indtrial(trial{j}); % and trials
Nk = length(c);
for k = 1:Nk
Y = A{i,m}*D{i}(Ic{i,m},It{i},c(k));
MY{m} = MY{m} + Y;
YY = YY + Y'*Y;
N = N + Nb(i);
end
end
% Apply any Hanning and filtering
%------------------------------------------------------------------
YY = W'*YY*W;
YY = T'*YY*T;
% Scale data (to remove subject and modality scaling differences)
%------------------------------------------------------------------
scale(i,m) = sign(trace(MY{m}'*(UL{m}*UL{1}')*MY{1}));
scale(i,m) = scale(i,m)/sqrt(trace(YY)/(Nm(m)*N));
YTY = YTY + YY*(scale(i,m)^2);
end
% temporal projector (at most Nmax modes) S = T*V
%======================================================================
if STAT % Stationarity assumption
S{i} = T; % temporal modes
Nr(i) = size(T,2); % number of temporal modes
VE(i) = 1; % variance explained
else
[U E] = spm_svd(YTY,exp(-8)); % get temporal modes
E = diag(E)/trace(YTY); % normalise variance
Nr(i) = min(length(E),Nmax); % number of temporal modes
S{i} = T*U(:,1:Nr(i)); % temporal modes
VE(i) = sum(E(1:Nr(i))); % variance explained
end
fprintf('Using %i temporal modes for subject %i, ',Nr(i),i)
fprintf('accounting for %0.2f percent average variance\n',full(100*VE(i)))
% whitening
%----------------------------------------------------------------------
Vq{i} = S{i}*inv(S{i}'*qV{i}*S{i})*S{i}'; % temporal precision
% get spatial covariance (Y*Y') for Gaussian process model.
%======================================================================
% loop over Nt trial types
%----------------------------------------------------------------------
UYYU{i} = 0;
for j = 1:Nt(i)
UY{i,j} = sparse(0);
c = D{i}.indtrial(trial{j});
Nk = length(c);
% loop over epochs
%------------------------------------------------------------------
for k = 1:Nk
% stack (scaled aligned data) over modalities
%--------------------------------------------------------------
for m = 1:Nmod
Y = D{i}(Ic{i,m},It{i},c(k))*S{i};
MY{m} = A{i,m}*Y*scale(i,m)/Nk;
end
% accumulate first & second-order responses
%--------------------------------------------------------------
Nn(i) = Nn(i) + Nr(i); % number of samples
Y = spm_cat(MY); % contribution to ERP
YY = Y*Y'; % and covariance
% accumulate statistics (subject-specific)
%--------------------------------------------------------------
UY{i,j} = UY{i,j} + Y; % condition-specific ERP
UYYU{i} = UYYU{i} + YY; % subject-specific covariance
% and pool for optimisation of spatial priors over subjects
%--------------------------------------------------------------
AY{end + 1} = Y; % pooled response for MVB
AYYA = AYYA + YY; % pooled response for ReML
end
end
end
% and concatenate for optimisation of spatial priors over subjects
%--------------------------------------------------------------------------
AY = spm_cat(AY); % pooled response for MVB
UL = spm_cat(UL); % pooled lead fields
% generate sensor error components (Qe)
%==========================================================================
AQ{1} = sparse(0);
for m = 1:Nmod
Qe{m} = sparse(0);
end
% assuming equal noise over subjects (Qe{m}) and modalities AQ
%--------------------------------------------------------------------------
N = cell(Nmod,Nmod);
for i = 1:Nl
for m = 1:Nmod
N{m,m} = sparse(Nm(m),Nm(m));
end
for m = 1:Nmod
Q = N;
AQeA = A{i,m}*QE{i,m}*A{i,m}';
Q{m,m} = AQeA/(trace(AQeA)/Nm(m));
Q = spm_cat(Q)/Nl;
Qe{m} = Qe{m} + Q;
AQ{1} = AQ{1} + Q;
end
end
%==========================================================================
% Step 1: Optimise spatial priors over subjects
%==========================================================================
% create source components (Qp)
%==========================================================================
switch(type)
case {'MSP','GS','ARD'}
% create MSP spatial basis set in source space
%------------------------------------------------------------------
Qp = {};
LQpL = {};
Ip = ceil([1:Np]*Ns/Np);
for i = 1:Np
% left hemisphere
%--------------------------------------------------------------
q = QG(:,Ip(i));
Qp{end + 1}.q = q;
LQpL{end + 1}.q = UL*q;
% right hemisphere
%--------------------------------------------------------------
[d j] = min(sum([vert(:,1) + vert(Ip(i),1), ...
vert(:,2) - vert(Ip(i),2), ...
vert(:,3) - vert(Ip(i),3)].^2,2));
q = QG(:,j);
Qp{end + 1}.q = q;
LQpL{end + 1}.q = UL*q;
% bilateral
%--------------------------------------------------------------
q = QG(:,Ip(i)) + QG(:,j);
Qp{end + 1}.q = q;
LQpL{end + 1}.q = UL*q;
end
case {'LOR','COH'}
% create minimum norm prior
%------------------------------------------------------------------
Qp{1} = speye(Ns,Ns);
LQpL{1} = UL*UL';
% add smoothness component in source space
%------------------------------------------------------------------
Qp{2} = QG;
LQpL{2} = UL*Qp{2}*UL';
case {'IID','MMN'}
% create minimum norm prior
%------------------------------------------------------------------
Qp{1} = speye(Ns,Ns);
LQpL{1} = UL*UL';
case {'EBB'}
% create beamforming prior. See:
% Source reconstruction accuracy of MEG and EEG Bayesian inversion approaches.
% Belardinelli P, Ortiz E, Barnes G, Noppeney U, Preissl H. PLoS One. 2012;7(12):e51985.
%------------------------------------------------------------------
InvCov = spm_inv(YY);
allsource = zeros(Ns,1);
Sourcepower = zeros(Ns,1);
for bk = 1:Ns
normpower = 1/(UL(:,bk)'*UL(:,bk));
Sourcepower(bk) = 1/(UL(:,bk)'*InvCov*UL(:,bk));
allsource(bk) = Sourcepower(bk)./normpower;
end
allsource = allsource/max(allsource); % Normalise
Qp{1} = diag(allsource);
LQpL{1} = UL*diag(allsource)*UL';
end
% augment with exogenous (e.g., fMRI) source priors in pQ
%==========================================================================
for i = 1:length(pQ)
switch(type)
case {'MSP','GS','ARD'}
%--------------------------------------------------------------
if isvector(pQ{i}) && length(pQ{i}) == Ns
Qp{end + 1}.q = pQ{i}(:);
LQpL{end + 1}.q = UL*Qp{end}.q;
else
errordlg('Using MSP(GS/ARD) please supply spatial priors as vectors')
return
end
case {'LOR','COH','IID','MMN'}
%--------------------------------------------------------------
if isvector(pQ{i}) && length(pQ{i}) == Ns
pQ{i} = pQ{i}(:);
Qp{end + 1} = sparse(diag(pQ{i}.^2));
LQpL{end + 1} = UL*Qp{end}*UL';
elseif size(pQ{i},1) == Ns && size(pQ{i},2) == Ns
Qp{end + 1} = pQ{i};
LQpL{end + 1} = UL*Qp{end}*UL';
else
errordlg('spatial priors are the wrong size')
return
end
end
end
if ~isempty(pQ)
fprintf('Using %d spatial source priors provided...\n',length(pQ));
end
% Inverse solution
%==========================================================================
QP = {};
LQP = {};
LQPL = {};
% Get source-level priors (using all subjects)
%--------------------------------------------------------------------------
switch(type)
case {'MSP','GS'}
% Greedy search over MSPs
%------------------------------------------------------------------
Np = length(Qp);
Q = sparse(Ns,Np);
for i = 1:Np
Q(:,i) = Qp{i}.q;
end
% Multivariate Bayes
%------------------------------------------------------------------
MVB = spm_mvb(AY,UL,[],Q,AQ,16);
% Accumulate empirical priors
%------------------------------------------------------------------
Qcp = Q*MVB.cp;
QP{end + 1} = sum(Qcp.*Q,2);
LQP{end + 1} = (UL*Qcp)*Q';
LQPL{end + 1} = LQP{end}*UL';
end
switch(type)
case {'MSP','ARD'}
% ReML - ARD
%------------------------------------------------------------------
[C,h] = spm_sp_reml(AYYA,[],[Qe LQpL],sum(Nn));
% Spatial priors (QP)
%------------------------------------------------------------------
Ne = length(Qe);
Np = length(Qp);
hp = h([1:Np] + Ne);
qp = sparse(0);
for i = 1:Np
if hp(i) > max(hp)/128;
qp = qp + hp(i)*Qp{i}.q*Qp{i}.q';
end
end
% Accumulate empirical priors
%------------------------------------------------------------------
QP{end + 1} = diag(qp);
LQP{end + 1} = UL*qp;
LQPL{end + 1} = LQP{end}*UL';
end
switch(type)
case {'IID','MMN','LOR','COH','EBB'}
% or ReML - ARD
%------------------------------------------------------------------
Q0 = exp(-2)*trace(AYYA)/sum(Nn)*AQ{1}/trace(AQ{1});
[C,h] = spm_reml_sc(AYYA,[],[Qe LQpL],sum(Nn),-4,16,Q0);
% Spatial priors (QP)
%------------------------------------------------------------------
Ne = length(Qe);
Np = length(Qp);
hp = h((1:Np) + Ne);
qp = sparse(0);
for i = 1:Np
qp = qp + hp(i)*Qp{i};
end
% Accumulate empirical priors
%------------------------------------------------------------------
QP{end + 1} = diag(qp);
LQP{end + 1} = UL*qp;
LQPL{end + 1} = LQP{end}*UL';
end
%==========================================================================
% Step 2: Re-estimate for each subject separately (fusing all modalities)
%==========================================================================
for i = 1:Nl
fprintf('Inverting subject %i\n',i)
% generate sensor component (Qe) per modality
%----------------------------------------------------------------------
AQ = 0;
Qe = {};
for m = 1:Nmod
N{m,m} = sparse(Nm(m),Nm(m));
end
for m = 1:Nmod
Q = N;
AQeA = A{i,m}*QE{i,m}*A{i,m}';
Q{m,m} = AQeA/(trace(AQeA)/Nm(m));
Qe{m} = spm_cat(Q);
AQ = AQ + Qe{m};
end
% using spatial priors from group analysis
%----------------------------------------------------------------------
Np = length(LQPL);
Ne = length(Qe);
Q = [Qe LQPL];
% re-do ReML (with informative hyperpriors)
%======================================================================
Q0 = exp(-2)*trace(UYYU{i})/Nn(i)*AQ/trace(AQ);
[Cy,h,Ph,F] = spm_reml_sc(UYYU{i},[],Q,Nn(i),-4,16,Q0);
% Data ID
%----------------------------------------------------------------------
ID = spm_data_id(AYYA);
% Covariance: sensor space - Ce and source space - L*Cp
%----------------------------------------------------------------------
Cp = sparse(0);
LCp = sparse(0);
hp = h(Ne + (1:Np));
for j = 1:Np
Cp = Cp + hp(j)*QP{j};
LCp = LCp + hp(j)*LQP{j};
end
% MAP estimates of instantaneous sources
%======================================================================
M = LCp'/Cy;
% conditional variance (leading diagonal)
% Cq = Cp - Cp*L'*iC*L*Cp;
%----------------------------------------------------------------------
Cq = Cp - sum(LCp.*M')';
% evaluate conditional expectation (of the sum over trials)
%----------------------------------------------------------------------
SSR = 0;
SST = 0;
J = {};
for j = 1:Nt(i)
% trial-type specific source reconstruction
%------------------------------------------------------------------
J{j} = M*UY{i,j};
% sum of squares
%------------------------------------------------------------------
SSR = SSR + sum(var((UY{i,j} - UL*J{j}),0,2));
SST = SST + sum(var( UY{i,j},0,2));
end
% accuracy; signal to noise (over sources)
%======================================================================
R2 = 100*(SST - SSR)/SST;
fprintf('Percent variance explained %.2f (%.2f)\n',full(R2),full(R2*VE(i)));
% Save results (for first modality)
%======================================================================
inverse.type = type; % inverse model
inverse.smooth = s; % smoothness (0 - 1)
inverse.xyz = xyz; % VOI (XYZ)
inverse.rad = rad; % VOI (radius)
inverse.scale = scale(i,:); % data scale-factor
inverse.M = M; % MAP projector (reduced)
inverse.J = J; % Conditional expectation
inverse.Y = UY(i,:); % ERP data (reduced)
inverse.L = UL; % Lead-field (reduced)
inverse.qC = Cq; % spatial covariance
inverse.qV = Vq{i}; % temporal correlations
inverse.T = S{i}; % temporal projector
inverse.U = A(i,:); % spatial projector
inverse.Is = Is; % Indices of active dipoles
inverse.It = It{i}; % Indices of time bins
inverse.Ic = Ic(i,:); % Indices of good channels
inverse.Nd = Nd; % number of dipoles
inverse.pst = pst{i}; % peristimulus time
inverse.dct = dct{i}; % frequency range
inverse.F = F; % log-evidence
inverse.ID = ID; % data ID
inverse.R2 = R2; % variance explained (reduced)
inverse.VE = R2*VE(i); % variance explained
inverse.woi = w{i}; % time-window inverted
inverse.modality = modalities; % modalities inverted
% save in struct
%----------------------------------------------------------------------
D{i}.inv{D{i}.val}.inverse = inverse;
D{i}.inv{D{i}.val}.method = 'Imaging';
% and delete old contrasts
%----------------------------------------------------------------------
try
D{i}.inv{D{i}.val} = rmfield(D{i}.inv{D{i}.val},'contrast');
end
% Display
%======================================================================
if ~spm('CmdLine'), spm_eeg_invert_display(D{i}); end
end
if length(D) == 1, D = D{1}; end
c4science · Help