bug in mTRFtransform.m

[nbara]

Hi Mick,

There might be a small bug in mTRFtransform.m:

Consider the following toy example:

% Toy dataset
fs = 8;
X = rand(8, 500) / 4 - 0.125;  % eeg
y = rand(1, 500) / 4 - 0.125;  % stim
tmin = -125;
tmax = 250;

% Add some here and there, make sure spikes are delayed in X relative to y
X(3:end, [102, 202, 302, 402]) = 3;
y(:, [100, 200, 300, 400]) = 4;
X(1, :) = 0;  % zero out channel 1

% Train a backward model and invert it
[g,t,c] = mTRFtrain(y',X',fs,-1,tmin, tmax,1e-3);
[gt,tt,ct] = mTRFtransform(y',X', g, fs, -1, tmin, tmax, c);

% Gymnastics to format the new model weights
gtt = zeros(size(gt, 2), 1, size(gt, 1));
gtt(:, 1, :) = gt';
recon = mTRFpredict(y', [], gtt, fs, 1, tmin, tmax, ct');

% Plot
figure(1); clf;
plot(y, 'k:'); hold on;
plot(recon(:,3) - mean(ct));

There's an offset between the stim and the reconstruction:

I think maybe this line should be changed from :
model_t = reshape(model_t(size(x,2)+1:end,:),size(x,2),length(tmin:tmax),size(y,2));
to :
model_t = reshape(model_t(1:end-size(x,2),:),size(x,2),length(tmin:tmax),size(y,2));

EDIT1: upon closer inspection of the code my solution is probably not correct either. Not sure where this comes from

[mickcrosse]

Hi Nicolas,

Thank you for flagging this issue with us. I made substantial changes to mTRFtransform after I read your comments but forgot to respond to you. There were indeed several issues: 1) it was not compatible with multivariate stimulus data and 2) was not outputting the transformed model with the same shape as a forward model (as I think you pointed out).

The reason that there is an offset between the stim and the reconstruction in your code is because you have not reconstructed the stimulus, you have reconstructed (or rather predicted) the EEG instead. Remember, you have used a transformed backward model, which is now functioning like a forward model, and you have convolved it with the stimulus y, giving a prediction of the EEG X. Thus the offset is the delay you created yourself between y and X. If instead you plot your reconstruction with the EEG X, you'll see they match up perfectly.

I have attached a modified version of your code as some of the functions have now changed in terms of their input/output arguments. I hope this helps and let me know if you've any questions.

Cheers,
Mick

% Toy dataset
stim = rand(500,1)/4 - 0.125;
resp = rand(500,8)/4 - 0.125;

% Params
fs = 8;
Dir = -1;
lambda = 1e-3;
tmin = -125;
tmax = 250;

% Add some spikes and make sure they are delayed in resp relative to stim
resp([102, 202, 302, 402],3:end) = 3;
stim([100, 200, 300, 400],:) = 4;
resp(:,1) = 0;  % zero out channel 1

% Train a backward model
bmodel = mTRFtrain(stim,resp,fs,Dir,tmin,tmax,lambda);

% Transform to a forward model
model = mTRFtransform(bmodel,resp);
model.b = repmat(bmodel.b,1,8);

% Predict EEG
pred = mTRFpredict(stim,resp,model);

% Plot EEG and prediction
figure(1); clf;
plot(resp(:,3),'k:'); hold on;
plot(pred(:,3));