trouble getting same residuals with BoostingResult and Convolve

[iranroman]

Hello eelbrain,

I've been trying to understand the boosting objects and operations. A sanity check I want to achieve is calculating the residual by hand by comparing the output of convolve and the "target" signal that boosting was optimized to model. However, I have been unsuccessful. Below is a minimal example to replicate my logic (based on one of the repo's examples).

# Original author: Christian Brodbeck <christianbrodbeck@nyu.edu>
# modifications by Iran R. Roman
# sphinx_gallery_thumbnail_number = 4
import os

from scipy.io import loadmat
import mne
from eelbrain import *
import numpy as np

# Load the mTRF speech dataset and convert data to NDVars
root = mne.datasets.mtrf.data_path()
speech_path = os.path.join(root, 'speech_data.mat')
mdata = loadmat(speech_path)

# Time axis
tstep = 1. / mdata['Fs'][0, 0]
n_times = mdata['envelope'].shape[0]
time = UTS(0, tstep, n_times)
# Load the EEG sensor coordinates (drop fiducials coordinates, which are stored
# after sensor 128)
sensor = Sensor.from_montage('biosemi128')[:128]
# Frequency dimension for the spectrogram
band = Scalar('frequency', range(16))
# Create variables
envelope = NDVar(mdata['envelope'][:, 0], (time,), name='envelope')
eeg = NDVar(mdata['EEG'], (time, sensor), name='EEG', info={'unit': 'µV'})
spectrogram = NDVar(mdata['spectrogram'], (time, band), name='spectrogram')
# Exclude a bad channel
eeg = eeg[sensor.index(exclude='A13')]

# model eeg by boosting the speech envelope
res = boosting(eeg, envelope, -0.100, 0.400, basis=0.100, partitions=4)

# convolve the filters with the envelope to obtain the model prediction
prediction = convolve(res.h_source, envelope)

# calculate the residual error "by hand"
res_error_hand = np.sum((prediction.get_data(['time','sensor'])-eeg.get_data(['time','sensor']))**2,axis=0)

# assert that res.residual and the residual obtained by hand are the same
assert np.allclose(res_error_hand, res.residual.get_data(['sensor']))

Maybe there's something fundamental I'm not understanding. Thanks!

[christianbrodbeck]

You’d want to use h_scaled to apply to the original data, not h_source.

[iranroman]

Hello @christianbrodbeck . Thanks a lot for your response. I made the change you suggested, but still cannot get the assertion to be true. What else could I be doing wrong?

# Author: Christian Brodbeck <christianbrodbeck@nyu.edu>
# sphinx_gallery_thumbnail_number = 4
import os

from scipy.io import loadmat
import mne
from eelbrain import *
import numpy as np

# Load the mTRF speech dataset and convert data to NDVars
root = mne.datasets.mtrf.data_path()
speech_path = os.path.join(root, 'speech_data.mat')
mdata = loadmat(speech_path)

# Time axis
tstep = 1. / mdata['Fs'][0, 0]
n_times = mdata['envelope'].shape[0]
time = UTS(0, tstep, n_times)
# Load the EEG sensor coordinates (drop fiducials coordinates, which are stored
# after sensor 128)
sensor = Sensor.from_montage('biosemi128')[:128]
# Frequency dimension for the spectrogram
band = Scalar('frequency', range(16))
# Create variables
envelope = NDVar(mdata['envelope'][:, 0], (time,), name='envelope')
eeg = NDVar(mdata['EEG'], (time, sensor), name='EEG', info={'unit': 'µV'})
spectrogram = NDVar(mdata['spectrogram'], (time, band), name='spectrogram')
# Exclude a bad channel
eeg = eeg[sensor.index(exclude='A13')]

# model eeg by boosting the speech envelope
res = boosting(eeg, envelope, -0.100, 0.400, basis=0.100, partitions=4)

# convolve the filters with the envelope to obtain the model prediction
prediction = convolve(res.h_scaled, envelope)

# calculate the residual error "by hand"
res_error_hand = np.sum((prediction.get_data(['time','sensor'])-eeg.get_data(['time','sensor']))**2,axis=0)

# assert that res.residual and the residual obtained by hand are the same
assert np.allclose(res_error_hand, res.residual.get_data(['sensor']))

[christianbrodbeck]

The residual (as all fit metrics) is calculated after normalization, so you'd have to take that normalization into account.

[iranroman]

Hello @christianbrodbeck , thanks for the reply.

What should be normalized and at what stage? Along which dimensions/features? Thanks for clarifying

[christianbrodbeck]

See scale_data boosting parameter. I’ll add this, will this be clear:

scale_data : bool | 'inplace'
    Scale ``y`` and ``x`` before boosting: subtract the mean and divide by
    the standard deviation (when ``error='l2'``) or the mean absolute
    value (when ``error='l1'``). Use ``'inplace'`` to save memory by scaling
    the original objects specified as ``y`` and ``x`` instead of making a 
    copy. Scale is stored in the :class:`BoostingResult: :attr:`.y_mean``,
    :attr:`.y_scale`, :attr:`.x_mean`, and :attr:`.x_scale` attributes.

[iranroman]

Thank you @christianbrodbeck . I think now I better understand what is going on. I assume that the whole point of getting both h and h_scaled is so that you can directly apply h_scaled to the original x. In contrast, the scale would have to be manually applied to y. I tried just that, but the numbers are still off. Maybe I'm still doing something wrong.

# Author: Christian Brodbeck <christianbrodbeck@nyu.edu>
# sphinx_gallery_thumbnail_number = 4
import os

from scipy.io import loadmat
import mne
from eelbrain import *
import numpy as np

# Load the mTRF speech dataset and convert data to NDVars
root = mne.datasets.mtrf.data_path()
speech_path = os.path.join(root, 'speech_data.mat')
mdata = loadmat(speech_path)

# Time axis
tstep = 1. / mdata['Fs'][0, 0]
n_times = mdata['envelope'].shape[0]
time = UTS(0, tstep, n_times)
# Load the EEG sensor coordinates (drop fiducials coordinates, which are stored
# after sensor 128)
sensor = Sensor.from_montage('biosemi128')[:128]
# Frequency dimension for the spectrogram
band = Scalar('frequency', range(16))
# Create variables
envelope = NDVar(mdata['envelope'][:, 0], (time,), name='envelope')
eeg = NDVar(mdata['EEG'], (time, sensor), name='EEG', info={'unit': 'µV'})
spectrogram = NDVar(mdata['spectrogram'], (time, band), name='spectrogram')
# Exclude a bad channel
eeg = eeg[sensor.index(exclude='A13')]

# model eeg by boosting the speech envelope
res = boosting(eeg, envelope, -0.100, 0.400, basis=0.100, partitions=4)

# convolve the filters with the envelope to obtain the model prediction
prediction = convolve(res.h_scaled, envelope)

# calculate the residual error "by hand"
res_error_hand = np.sum((prediction.get_data(['time','sensor'])-(eeg.get_data(['time','sensor'])-res.y_mean.get_data(['sensor']))/res.y_scale.get_data(['sensor']))**2,axis=0)

# assert that res.residual and the residual obtained by hand are the same
assert np.allclose(res_error_hand, res.residual.get_data(['sensor']))

[iranroman]

Also, given what we have discussed so far, I would have expected these to give the same result, but they don't.

# convolve the filters with the envelope to obtain the model prediction
prediction_hand = convolve(res.h, (envelope-res.x_mean)/res.x_scale)
prediction = convolve(res.h_scaled, envelope)

assert np.allclose(prediction_hand.get_data(['sensor','time']),prediction.get_data(['sensor','time'])*1000)

I feel like there's something very fundamental that I'm missing.

[proloyd]

Hi @iranroman,

Prediction

according to my understanding, the prediction should be:

prediction = convolve(res.h_scaled, (envelope - res.x_mean)) + res.y_mean

and that matches with hand calculated prediction:

prediction_hand = convolve(res.h, (envelope - res.x_mean) / res.x_scale) * res.y_scale + res.y_mean
assert np.allclose(prediction, prediction)

Residual

However, according to documentation, the residual should be:

error = ((eeg - res.y_mean) - convolve(res.h_scaled, (envelope - res.x_mean))) / res.y_scale
residual = (error ** 2).sum('time')

but this does not match with the res.residual. May be this residual is the residual from the data fitting phase (evaluated separately on the partitions and averaged) @christianbrodbeck?

[iranroman]

Hello, thanks you @proloyd for your response.

This line seems rather complicated to me:

prediction_hand = convolve(res.h, (envelope - res.x_mean) / res.x_scale) * res.y_scale + res.y_mean

Let me try to break this down:

convolve(res.h, (envelope - res.x_mean) / res.x_scale)

This expression convolves the unscaled TRF filters with the standardized envelope (?). I'm assuming here that res.x_scale is the envelope's standard deviation, but res.x_scale==0.0114 and np.std(np.array(envelope)))==0.0123, so not sure if I'm misunderstanding something there (please answer at least about this bit if you have time).

This is then multiplied by res.scale and res.y_mean are added in order to "unstandardize" the prediction and bring it back to the units of the original target (in this case the variable called eeg; by the way, here both np.mean(eeg.get_data(['time','sensor']),axis=0)==res.y_mean.get_data(['sensor'])
and np.std(eeg.get_data(['time','sensor']),axis=0)==res.y_scale.get_data(['sensor']) are True).

With all of this said, I would have expected a way to say something along the lines of convolve(res.h, envelope) to get the model's predictions in an "inference mode" kind of fashion. However, I can live with this.

P.D: I think in the second line of code you meant assert np.allclose(prediction, prediction_hand) instead of assert np.allclose(prediction, prediction)

[iranroman]

@christianbrodbeck , could you please provide me with an example for how to use BoostingResult in combination with convolve to apply it to a time-series signal and obtain a prediction that I can directly compare with EEG data? Thank you!

I'm trying to incorporate your boosting code as evaluation metric in the optimization routine of a model that predicts speech envelopes. Cheers!

[iranroman]

would it be using

convolve(res.h_scaled, (envelope - res.x_mean))*res.y_scale + res.y_mean

?

That's what would make sense to me, but without being able to replicate the residuals, I cannot corroborate the operations I'm trying to do "by hand". Thanks!

[proloyd]

would it be using

convolve(res.h_scaled, (envelope - res.x_mean))*res.y_scale + res.y_mean

?

That's what would make sense to me, but without being able to replicate the residuals, I cannot corroborate the operations I'm trying to do "by hand". Thanks!

What you want is:
prediction = convolve(res.h_scaled, (envelope - res.x_mean)) + res.y_mean

This aligns with the fact that res.h_scaled = res.h * res.y_scale / res.x_scale.

[christianbrodbeck]

How big is your difference? Could it be due to numerical precision?

For example, you're calculating it at a different scale than the actual code, and the SS is calculated using eelbrain._trf._boosting_opt.l2

For reference, the actual calculation is here:

Eelbrain/eelbrain/_trf/_boosting.py

Line 761 in e81c439

convolve(h[i_y], self.data.x, self.data.x_pads, self._i_start, segments, y_pred_i)

The L2 evaluator is here:

Eelbrain/eelbrain/_trf/_fit_metrics.py

Line 74 in e81c439

class L2(Evaluator):

[iranroman]

I tried prediction = convolve(res.h_scaled, (envelope - res.x_mean)) + res.y_mean and this is how the prediction compares visually to the original EEG. I could be wrong, but the plots look like they are off on the order or magnitude. Or does this look right to you @christianbrodbeck @proloyd ? Thanks again!

[christianbrodbeck]

Hard to say from this, boosting is pretty conservative. Notice that the visible features of the response are not in the prediction, so they are probably noise not signal.