Stephanie Haro: yeah. Do you want to start?
Ole: Yeah. So I'll start
Ole: sharing my screen.
Ole: Oh, you have to allow me.
Stephanie Haro: Let me give you preferences.
Stephanie Haro: The
Ole: you're a co-host now.
Ole: Okay. Can you see the notebook? Now.
Stephanie Haro: Yep.
Ole: okay, cool. So hello, everyone again.
Ole: So
Ole: II would like to make this. Wait. I first. I can see the chat like that I first have to go back here. So II wanna make it as interactive as possible. So whenever you at any point feel like something's unclear, or you have questions, or whatever. Just feel free to interrupt me. I want to do like a little
Ole: mini questionnaire, just to know how much of experience you have with the topics being covered.
Ole: So could you in the chat like type one, if you already did any kind of analysis, using a Trf or like a Trf model, and doesn't matter if it's in Python or Matlab and 2, if you never did any trf analyses yourself.
Ole: Okay, so kind of mix here. Most most people.
Ole: most people already have done some kind of trf analysis. And then again, the same thing. Maybe type one. If you feel somewhat confident in using python and 2, if you feel like you have very little experience using python.
Ole: Okay? Seems more more. Choose for that answer, I say, so. So yeah, just so. I know that I don't like spend too much time on extremely basic things. And but again, if feel free to to interrupt me at any time.
Ole: and so I'll get started. Now
Ole: I'm all wait. I have to share
Ole: alright.
Ole: So
As I already said in the talk, we'll be working with this data set of the Old Man and the sea.
Ole: I recently uploaded that to Openuro.
Ole: So if you've never used open euro, it's really cool. It's this online platform where people share data sets and the ho, all the data sets there are in the bids format, which is a standardized format for for storing data.
Ole: And the whole thing interacts with a software that's called data led
Ole: and who hasn't done it yet. Please execute the code here in all those sections up to
Ole: up to section number 2, because that will start the downloading of the of the data which takes a bit.
Ole: So so data led is a really cool piece of software for kind of managing data sets and analyses. III don't wanna dive into the details in this talk, because it's a a very deep
Ole: piece of software, and it takes a while to understand. But it's II encourage everyone to have a look at it. I think it has my made my life as a researcher, a lot better. II never now do any research projects where I don't use data, let so definitely take a look at their their online documentation and see if it is for you. For now we just needed to download the data from from open euro.
Ole: So here in this first cell, we're installing some software which is we're installing first data led.
Ole: And then we're using data led to install the data set from open euro. So this is the address of the data set. and you can see like, if I execute that section, it
Ole: runs really quickly. It's only a few like 3 s, and of course, in 3 s you're not downloading a whole data set. But instead, what data led does is it? It downloads the. It loads the whole structure of the data set. So
Ole: we see here. They're like all the all the folders.
Ole: like 2019 subjects in total. A folder with stimuli and everything. So it it downloads, though, and if you look into those folders, you will actually see all the files, for all the subjects are there. But those are not the actual files, or those are just symbolic links. And what's cool about this is that it's a way where you can very quickly get the whole data set and just work on parts of it.
Ole: And then push back your changes to your original repository. And that way you could still like kind of preserve the structure. Why? And you don't have to. Don't have to download the whole data set. This is really useful if you work with stuff like stuff from, for example, Human Connectome Project, where, of course, you you can't download the whole data set on your on your on your notebook drive
Stephanie Haro: all. Do you mind increasing your font size.
Stephanie Haro: We have that comment from the chat.
Stephanie Haro: Is that is that better
Stephanie Haro: on the chat?
Stephanie Haro: Let's go with yeah. Great.
Ole: great. Yeah. II can't see the chat. So if you you can just and tell me, or or you tell me, yeah, thanks.
Ole: And so if you want, if and if you're not familiar with how the bids data structure looks like you can just take a moment to
Ole: to look at the content of the different folders. So you can just use this command. Ls.
Ole: With some folders. So I've I've called the folder here all mats. So for old man and the C.
Ole: And you can
Ole: list out contents by, for example, saying. and something like
Ole: homeats. Subjects.
Ole: Okay, you're even getting the suggestions. And you can see each fold. Each subject folder has only a single sub folder called eg.
Ole: Because bits also
Ole: includes other modalities. It's like a a general neuroscience data type.
Ole: And then, if we look in that. Eg, folder.
Ole: we can see there's a bunch of files. So this is the brain vision data format meaning that for each file there is a, there's a Eg file. And then there's a a Json file. And then these like.
Ole: Vhdr and vm, rk, so this brain vision data set like splits the data across multiple files, and then for each
Ole: for each file. There is also a file called Channels, which lists the bad channels in that segment.
Ole: and the data is divided into runs. So each one run here means like one just continuous presentation of data.
Ole: And then here, which by now hopefully has finished for everyone. We're actually loading the data. So this is with data led, you can say there's a command called, Get
Ole: and the. And you can call this. Get command
Ole: on on any
Ole: folder.
Ole: and then data let will check. If the the contents in that folder are present on your machine, and if not download them. So here I select the subject, or 3, you can use a different subject. So you can just change this number, and you will be able to do all the analyses with a data of a different subject. If you, if you think that might be more interesting if you're feeling a little adventurous. 2.
Ole: Of course. Then it might be that you pick one of the subjects that has
Ole: which have really noisy responses. And you won't be able to see very nice results. So it's it's a bit of a dangerous choice if you just select at random.
Ole: but feel free to do that if you want.
Ole: And then finally, we have to install some python dependencies where, of course, use Mtrf, that that's the package we're interested in. I also download M and E, which we have to use for loading and Pre processing the Eg data. And then another small package called S lab, which is another project I'm involved when
Ole: and which is a toolbox for psycho acoustic acoustics and sound processing, which here I'm just using to to load and pre-process audio
Ole: but if you do psycho acoustics you might wanna wanna check that out. II put a link in here.
Ole: Okay, so did that work for everyone. Anyone who has had trouble downloading the data.
Ole: No, that's that's great.
Ole: So we'll start with a a simple, univariate model of the acoustic envelope. So probably, if you're all familiar with Trf analyses, you know that the envelope is probably the most widely used or or most studied feature when it comes to speech processing. There is also a bunch of reviews on on
Ole: on envelope tracking studies. I I've linked one here, for example. And so first we have to compute the envelope. So we oh, here, I think I didn't mention that up here
Ole: with the
Ole: with the get command. We also downloaded the contact, content content of the stimuli folder which contains the the waf files of the presented audio, and also word and phony annotations which will not use here. But if you're interested in doing something with that, they're also in there.
Ole: And so first, we have to compute the envelope
Ole: and the envelope is obtained by that's maybe for start the code, because it it takes quite a while to run. This is, I think, the the
Ole: bit of the code that is takes the longest to run of. Of all the stuff in here. The envelope is obtained by
Ole: by computing the the Hilbert transform of the data
Ole: and then so the which is gives you the analytic signal.
Ole: And then you down sample that to 64 Hertz, because usually we're interested in the low, frequent part of the eg. Response. And so there's no reason for us to keep to keep all the data to resolve high frequencies. So down sampling everything to 64, Hertz will will make it a lot quicker to compute.
Ole: One thing that's really important here is to to down sample to 64 hertz after computing the envelope not before. Because if you down sample to 64 hertz
Ole: immediately, you'll lose most of the information of the speech signal, because most stuff happens above 64 hertz, and then you'll still get out an envelope. That kind of looks like an envelope should look like, but it won't carry any information. And actually, that recently happened to me where I had some analysis running and the
Ole: just, the all the predictions were complete garbage, and it took me a whole day due to notice that I did be resampling before computing the envelope. And and that was actually the issue.
Ole: So I am down sampling here before and before doing the envelope, but only to
Ole: 1,400 1,470 hertz. So this is a and just to like speed up the computation. So this is still enough high of a sampling rate to get most of the interesting stuff in the speech signal.
Ole: and we can see, okay, it's it's still running on.
Ole: so maybe we can already look at the plot. You you should be able to see the plot here? Right?
Ole: Or is it just for me?
Stephanie Haro: Are you also seeing the plots in your own notebooks? Yeah.
Ole: okay,
Ole: So
Ole: here. I just pick like a short segment, and and plot the speech and the envelope, and we can.
Ole: And you can see that the envelope is kind of just
Ole: similar, as you're like drawing an outline along the speech.
Ole: and we'll use that for like now, the first part of the or for actually, for the most part of this tutorial. We'll use the envelope as our predictor
Ole: does it look similar for everyone else? You also have, like
Ole: half of the files processed.
Stephanie Haro: Could you do a little survey with like a a one? If you're
Stephanie Haro: processing the the stimulus envelope
Stephanie Haro: or or 2? Or do you want to coordinate that
Stephanie Haro: survey
Stephanie Haro: like, are some people still downloading the Eeg.
Stephanie Haro: yeah, that's useful.
Ole: So is that all? Yeah. Envelope? Okay? Similar
Ole: similar to to what I'm experiencing, and maybe we should have also ran that earlier. But II promise the remainder will not take so long. And then we after that we have another another part that that takes a bit maybe I I'll go over this now. So then,
Ole: After after we got the envelope. We'll load the eg. Data and section 2.2
Ole: erez agmoni, and this uses M. And E. Python. Oh, and maybe a one, if you're familiar with M. And E. Python, and a 2 if you're not 150.
Ole: Okay. most people haven't used them in E. Python.
Ole: Yeah. So M. And E is kind of the the most widely used python toolbox for working with. Eg, data. It's it's very nice. I like it. A lot has great online documentation.
Ole: so here, what we're doing is we're listing all the we're listing all the files. So I'm not going to. I'm I'm if no one asks about it. I'm not going to explain any like one particular
Ole: python, command and detail, but rather just kind of explain what they're supposed to achieve. If you have like particular questions, feel free to ask
Ole: mit Ctl, and so we're here creating a electrode montage because our our eg. Data does not contain information about channel locations, but will will need that information for the pre-processing. So we just create the standard montage for the 128 Channel Biosmicap, which is the one that was used to record the data.
Ole: And then here we'll just list out all the files
Ole: and and sort them
Ole: in this line here, and then I'll print them out to see that they're in the right order, because that's sometimes can be tricky, and when you mess up the order of your files the the results can be pretty bad. But here it seems, okay. So we have like.
Ole: run one run 2. That's why it's important if you number them, that you
Ole: don't number them. Just 1, 2, but oh, 1, 0 2, because,
Ole: otherwise the the the order will be messed up.
Ole: and then we go through all those recordings
Ole: and and channel files. So the channel files just contain the listings of bad channels. We. We load the the data using these M and E functions?
Ole: we interpolate the the bad channels
Ole: then we filter between one and 15 hertz, which is also kind of the standard window people use for for envelope tracking.
Ole: and then we resample the Eg. Data to 64 hertz, just as we did with the envelope.
Ole: And then we then we re-reference the data to the global average. And then finally, we'll get the We. So this is like a an instance of the M. And E. Python, raw data class.
Ole: But we don't need this class. We just wanna get the numpy array inside the class that contains the data. So we use this get data method to just get the numpy array that's inside this raw object.
Ole: And then we have to transpose it. Meaning. We have to flip it around because m and E has the channels as their first dimension, and and time points as their second dimension. But most other stuff, and and the Mtrf.
Ole: Pi toolbox uses it the other way around, meaning the first dimension or the in python. It's Z starts counting at 0. So the 0 dimension should be time and the dimension one
Ole: should be should be channels. So we have to transpose it. And then it's just appended to this list. So right here in at the top, I'm creating an empty list. And then I'm I'm appending
Ole: at each
Ole: the raw file, the the data matrix of each recording to that list. So let's check.
Stephanie Haro: I have a basic question. And then, camila, if you wanna unmute yourself after asking a question. You could or I can ask your question. Just let me know. My basic question is, how difficult is it to
Stephanie Haro: take Eeg data that you parse out in Matlab and saved as a dot mat, and then set it up so it can be read. But set it up so it's an MNE.
Stephanie Haro: Data structure format.
That shouldn't be too hard. There is
Ole: so you can load Matt files just in python. There is the library. Scipy has a function that's called a load, Matt. So I think that's that should be
Stephanie Haro: like is the M and E data structure
Ole: analogous to a Matlab structure. So
Stephanie Haro: like, is it easiest to save a structure
Stephanie Haro: as a dot mat and then load it in as sci-fi you like. Yeah, I mean, if you, if you're just interested in using Mtrf. Pi, you don't need to use M. And E at all.
Ole: I mean Emini
Ole: with M and E, you can use all their like nice pre processing functions like filtering, and so on. But if you don't want to use that, you don't have to use M. And E, because in the end, right, we're anyway throwing away the M. And E object, and just keep the numpy array so and I think if you do like loadmat, you're you're also getting out an umpire ray, so you can also
Ole: just use that and then use some other functions for filtering and and not use M and E functions, but it also in M and E, you can just create a new
Ole: raw data object from scratch where you
Ole: yeah, here is like a tutorial on that where you it's a bit of a longer process, because you have to add, like manually add all this Meta information like the names of channels, and so on, that probably are not included in the in the mad file.
Ole: Okay.
Ole: so it might be easier to just skip the the M. And E. Step.
Stephanie Haro: Camila, would you like to unmute yourself to ask your question?
Camila Zugarramurdi: Yeah, thanks. My question was regarding the reference because I saw you were using the global average. And I was wondering whether, for computing Trfs in particular, is there any recommendation regarding which reference to use, or it doesn't really matter.
Ole: Yeah, I think I think the generally in all Eeg stuff. The question of reference is like really a tough one, and there's
Ole: it's a very deep literature, and I don't know it has been debated for since eg. Is around, what's the right reference? And people have different opinions on it.
Ole: I like to use the average as kind of a default. If you don't have a compelling reason not to, because it's it's robust meaning. It doesn't really like a single bad channel doesn't have a a big effect. For example, people often use the link mastoids.
Ole: which it is fine, but a downside of that is, if, then, one of your master channels is bad, this these artif, and has some artifacts in them that they will project into all the other channels. And this is less of a problem with the average. And then the average is also kind of
Ole: more fair to like different kinds of signals. With different topographies. If you're I mean, if you know that you're only looking for auditory stuff, it's fine using a mastoid. But II think the average is like a reasonable default. If you don't have any any strong opinions about any particular reference.
Thank you.
Ole: okay. So the envelope has computed. So now we can run those lines here.
Ole: this will also take a moment, but not as long as the the envelope.
Ole: And so from this we'll get
Ole: one list of responses, and that this should contain 20 elements, one for each run.
Ole: because there are 20 runs for each subject.
Ole: and each and each array in this in this list should have the shape. Samples by channels and channels are 128, and samples varies a bit between the segments.
Ole: And then, once we have that, we have to do one more thing.
Ole: which is, we have to crop the data to the same length, because usually the recordings are a bit longer than the stimulus. So after the the speech stops, the recording goes on for a bit. And so we we have to go through there and and crop them. I've written it here in the way that
Ole: you crop you check what whatever is is longer the stimulus and the response, and then crop the other one
Ole: to the, to the duration of the shorter one. It's for almost all subjects the responses is longer than the stimulus. There's like one or one subject
Ole: in the data set where the where the response is like a few samples shorter. So it seems like that recording will stop prematurely or something like that. Which is why I included that that statement.
Ole: And then we normalize the data by z-scoring, which means we subtract the mean and divide by the standard deviation. And that means that then, if we would plot the data, the Y
Ole: the Y's X's would tell us for a particular data point, how many standard deviations away from the mean, that is. And this is just a way to like a standard way for for linear modeling. And it's especially important if you have multiple features
Ole: that are on a different scale. Because you you don't, you wanna be able to compare the weights across features. But this won't be possible if they have, like, a widely different scale, because if you have like one feature with much bigger absolute values. The model will assign smaller weights because they
Ole: because
Ole: you have base, you have basically have to multiply the feature by the weight to get the to get the effect on the response, and if they, if you have like absolutely large features, you have to like, make the wait smaller to account for that. And so to to be basically have like a fair comparison of all features you usecore them here. In this case it it probably wouldn't matter, because we're just using a single feature. But just for sakes of like keeping the standard workflow.
Ole: I included the as he's scoring here.
Ole: so it seems the eg. Pre-processing has run. So we can do the Z scoring. Okay.
Ole: so now, to the top, we, we have a question is this at the time for a question. So with respect to the stimulus in the trials, is this like one single track of all the speech someone heard, or is it cut into different trials, and how long the trials are on?
Ole: Yeah, so the it is into, I mean, they're called runs here right? If you look at the list. It's like they're they're the files are called run one, run, 2, run 3. And II mean, II didn't record this data, but as far as I know, one run is like one continuous presentation of data. So you hear
Ole: you hear one run the the speech for one run. Then there's a short break, I think. People answered some comprehension questions, and then the next one is played, and so there, I think they're roughly like 2 min each. They vary a bit.
Ole: But yeah, they're like organized in in trials.
Ole: Does it answer your question?
Camila Zugarramurdi: Yes, thanks.
Ole: Alright.
Ole: Okay. So now,
Ole: to the time lags. So when we when we model the relationship between stimulus response.
Ole: we have to account for the fact that this relate that this, the response is not instantaneous meaning it takes. If you hear the stimulus, it takes a certain time until the stuff reaches the cortex, and we see something.
Ole: and to account for this, we compute the the Trf computes regression across multiple time lags.
Ole: And we put the time lab, we do this with a single match matrix by sort of stacking time lag copies of the stimulus.
Ole: So if you run so here, I'm using this from Mtrf dot matrices, I'm using this function lag matrix to generate this time. Lag matrix X. This is something. So the next few steps is stuff you you're usually not doing when you do an actual trf analysis. But I wanna kind of walk you through the sort of nuts and bolts of the the method. So you have, like a solid understanding of what's what's happening under the hood.
Ole: So if we execute this, we'll get a plot here. Where in the top plot you're seeing the stimulus feature, the envelope at multiple time legs here the blue one at 0, and then the the green one is something like 500 ms
Ole: and
Ole: and you can see in the beginning there is this
Ole: for the yellow, or orange, and green line. There is this step in the beginning here.
Ole: and and this happens because of the because of
Ole: because the time lagged features are padded with zeros. So they're like all of the same length, and they're on the left side. They're padded with zeros, and because we but because we z scored.
Ole: and the the baseline of our envelope is below 0. That is because before Z. Scoring the baseline of the envelope is 0, and then, if you subtract the mean, it, it gets shifted a bit, and so if you, if you then pad it with zeros. You're getting this step in the beginning, and this can cause sometimes artifacts at the edges of trfs. So it's good to know where where those come from.
Ole: And then here we see the neural response in the same time window, and you could see it's it's pretty noisy, and there's not really any
Ole: any relationship you can just see by by eyeballing it. So we have to. We have to use statistics to to figure out if there's any relationship.
Ole: Now, I have a little exercise, because I want to keep this interactive. And I'll just talk for 2 h.
Ole: And that is so here I plot the first 500 samples of the signal. and what I would like you to do is try to plot the last 500 samples of the of the time left envelopes.
Ole: So I'll give you a few moments to figure out how you would do that. And then, continue.
Stephanie Haro: I, subconscious is so Matlab trained that I googled Matlab Index last elements of vector, without even noticing.
Ole: Okay, I think I don't wanna wait too long here because I don't wanna just use up all our time for for array. Indexing and python has this really cool thing where you can, where you can index arrays from from the back
Ole: by putting a minus in front of the index. So I can say, Here X is the stimulus, the timelike matrix. So I can say, X, minus
Ole: 500, and then just leave the other
Ole: are here, empty, and leaving an empty just says, Take everything. So this says, take plot, E, take everything from 500 samples before the end, and then up to the end.
Ole: And then, if we run that.
Ole: we'll see now, here, this is the the last samples, and you can see that here now there's no 0 padding at this end, but rather the time. Lag copies are are cropped in the end.
Ole: and this is because there's just a The.
Ole: We can't quantify the the response to the envelope
Ole: at a later point than the the last sample of the unleagged envelope. Right? It's just the the question, what is the
Ole: what's the neural response that happens?
Ole: 300 ms after this last point of the blue line. Here is just something we we can't answer because we don't have any brain recording of that time. And so you're just at the other end. It's just cut off.
Ole: And so this is the the timeline matrix. And I just wanted to take a moment to to show you how that looks like
Ole: so now, once we have that time, like Matrix X, we can compute the ordinary least square regression which has the closed form solution that's given here, which is the X transposed, multiplied by X
Ole: and this is called the the auto covariance of the stimulus
Ole: it kind of tells you how strongly the stimulus is correlated with itself. and then this is inverted, shown here by the minus one. and this is then multiplied by the covariance of the stimulus and the response.
Ole: And the reason that we use this inverse auto covariance of the stimulus is because we wanna account for the auto covariance like we. We can't just consider the covariance of the stimulus and the response, because this is affected by the by the bias by the fact that the stimulus is correlated with itself.
Ole: And so the result of that this matrix W those are the weights of the trf, so this is a matrix that's shaped time lags by channels.
Ole: Or if you have multiple features, it would actually be a threed matrix that is, features by lags by channels. But since we only have a single one, it's just lags by channels.
Ole: and those are regression coefficients. So they they indicate how a unit change in the stimulus feature is expected to affect the neural response over time.
Ole: And first, let's take a look at the covariance matrices.
Ole: So here on the left you see the auto covariance of the stimulus.
This first term here is
Ole: this is for the for the a bias term in the
Ole: and the regression.
Ole: So the biases in linear regression you have like a plus x times B, one
Ole: plus x times b, 2, and so on. So the the the a is called the is the intercept or bias term and then the rest is just your your time lags.
Ole: and you can see there's a a diagonal where all the values are one. And this is just because, of course, every time like is correlated with itself, is identical to itself. And so the correlation is one, and then the more you go away from that the lower the covariance value gets. So this looks kind of like a ridge. If you would look at. Look at it from the top and imagine it. It's like the peaks along the diagonal, and then it falls off.
Ole: and the fact that it's not 0 off. The diagonal is indicative of the autocorrelation. And then, if we look to the right, we see the the covariance of
Ole: stimulus and response, and you can see there's some, and here it's time legs by channels, and you can see there's
Ole: difference across channels. But it looks very smeared across the time domain. And that is because.
Ole: due to the autocorrelation.
Ole: we can't infer like any temporal relationship from this. Just from this covariance matrix alone.
Ole: And this is why we're multiplying it with the inverse of that. And so now we have the matrices. And we can. We can run the linear regression.
Ole: And if we run that code, we get something really horrible like this kind of spiky looking, weird thing.
Ole: And so why does it look like that? The the reason we're getting something like this is because, all our
Ole: all our time lag the our different predictors, which here are like copies of the same signal at at different time lags very similar. And so that means that there are a lot of a lot of possible
Ole: solutions to the to the regression problem, meaning that you could. There, there would be a lot of different distributions of weights across the different time lags that would produce almost equally good solutions. Right? It. It probably doesn't matter if you use, like
Ole: one or the other of 2 neighboring time lags. And so you there's a lot of solutions to the to the regression fit that are almost equally good.
Ole: Can you explain what we're seeing on the y-axis? So the color is, each line is a different lug.
Ole: Oh, no. E, yeah. Each line is one channel. I'm sorry. Very good question. Yeah. Each line is 1 one eg. Channel that's being predicted, and the the amplitude shows you the weight of the Trf, so
Ole: that means how strongly the model thinks that a change in the envelope will affect the response at a given channel at a certain point in time
Ole: and
Ole: And so because they're like, there's many different solutions that are equally good. The the actual, like weight is more strongly determined by the noise than by the brain signal, which is why we see like this spiky looking thing. And this is a problem that's known as multi-colinarity.
Ole: If you wanna understand it more deeply. you can read this blog post. This is a very, very cool blog post that illustrates the problem.
Ole: But so basically, it's, there's a large space of possible solutions. And then the Trf, just kind of jumps around in that space. And but so to solve the problem of of colonelity.
Ole: we can use
Ole: regularization which will here use Ridge regularization, which is sort of the standard.
Ole: The standard way of of regular rising. T rf.
Ole: and this simply adds a small
Ole: or adds a constant to the to the auto covariance, which is the so
Ole: here, here have this factor, Lambda, and an identity matrix, and that is being added to the auto covariance before it's being inverted. And this has the the regularization has the effect that it that it penalizes large model weights. And so
Ole: and so this
Ole: this will like reduce the spikiness because it will it will dis encourage the Trf. To to take on those those extreme spikes. And so here I use the regularization value of like
Ole: 10 to the power of 6.
Ole: And now it looks very different. Now we're getting like much, much smoother curves. It still doesn't look very informative. But it is also because right here we're only using a I think
Ole: here you can see it.
Ole: I'm only using here one element of the stimulus list. So I'm only using a state single, a single trial of stimulus and response. So so it's not surprising that we get like a noisy
Ole: response. So
Ole: Now you can go ahead and try out a few different regularization values. And then just observe how the how the shape of the
Ole: functions changes as a function of of regularization.
Ole: I'll give you a moment for that. Yeah.
Camila Zugarramurdi: earlier in your talk you said that with higher lambda, you would get better solutions. So why, why, why wouldn't we always use like a high lambda? Well, it it depends. You're not always getting better solutions. There is.
Ole: There is, I mean, depending on on your on your data, and how much data you have and how noisy it is. The, there might be a different optimal value of lambda.
Ole: So that's a very good question. I'll I would
Ole: withhold the answer for a bit, because a later section is going to to address the question of how to choose the optimal lambda. And
Ole: for now we can just like play around a bit with it and see how the function changes.
Ole: Hmm.
Ole: and, for example.
Ole: I can set it to something much lower, and then.
Ole: when I set it to something low, it it again looks more like the original unregularized T. Rf.
Ole: and then when I when I
Ole: This was the value I had from the start, and then, when I set it to something much higher.
Ole: like 6 orders of magnitude higher. You can see it's
Ole: not really changing very much. So it seems like at some at some at some point.
Ole: It's the the regularization doesn't affect the shape of the Trf
Ole: too strongly.
Ole: The only thing that changes, if you look here is the scale of the weights. So here you can see right here. Now, we're at a really low scale, meaning the absolute values of the of the Trf are very small versus. When when I have this
Ole: this value here.
Ole: we get like higher, absolute values.
Ole: so this is the so with regularization, we're able to like, overcome this problem of multi-col linearity.
Ole: And so this was just supposed to be like a little overview of the the steps that happen
Ole: under the hood when you when you call, when you call the Ntrf. Train function. But so now that we've
Ole: we took a look at that.
Ole: we can just start using actually, the Mtrf pi functions, and make our lives a little easier. So instead of doing all that stuff that II did above.
Ole: we can just create a an instance of the Trf class.
Ole: and just I call it a forward Trf here.
Ole: and then we can just use the train function to train it on on our stimulus and response. Here, now we are using the whole data, not just a single trial. So we'll get out a nicer, a nicer trf, and then we can. We can plot the result.
Ole: And we see now we can get a like somewhat nicer looking trf.
Ole: and again, here each line is one channel, and you can interpret the the amplitude here as the expected change in the neural signal
Ole: following a unit change in the predictor. So if the if the if, the, for example, if we pick out that
Ole: that channel over here.
Ole: we're expecting, based on this model, that
Ole: if the envelope changes from 0 to one.
Ole: 200 ms later. We're expecting to see a change of
Ole: 0 point 2 at this channel.
Ole: although it. That is not entirely true. That that would be true if it's unregularized. But, as I shown you before, right? If I increase the regularization.
Ole: the scale changes so I can't
Ole: it. Once I start introducing regularization. I can't it. I can't
Ole: interpret the amplitude of the weights in absolute terms. So II can't. So that actually, my interpretation now is is not that I just gave you is not
Ole: exactly valid, because by regular rising we, we change the scale. However, I can still relatively compare the weights.
Ole: the channels. So what you can say is where, if following, if the the envelope changes from 0 to one.
Ole: we're expecting a a change in the neural signal after 200 ms, and that we're expecting that this changes like roughly twice as big at that channel than it is at that channel, and we're expecting that is positive at this channel, and we're expecting that it is negative at this channel. So do you. You can compare the the magnitude of weights
Ole: across. it across different channels and also across different trfs. If they use the same value for for regularization.
Camila Zugarramurdi: 2 questions, one. The the training function uses the whole dataset, or leaves some of the data out, and the second one is.
Camila Zugarramurdi: and the demine the the starting time.
Ole: Should it be like negative, so that you can see that there are no effects prior to stimuli.
Ole: That's actually what I have here in this exercise. So
Ole: yeah, feel free to to play around a bit with the team in and the team X
Ole: and
Ole: change it a bit, and then maybe think about what would be the appropriate and time window for this model.
Ole: Oh, sorry. I need to get some water.
Ole: Yeah, sorry. I'm still recovering from Covid, so there will probably be more coughing in this talk.
Stephanie Haro: Feel free to. If you need, like, a 5 min break, let us know when there's a good spot for everyone to take a bathroom break or have a 5 min break.
Ole: So
Ole: if we increase the team in T. Max Time window.
Ole: we can see that
Ole: there seems to be stuff happening before 0, although that is very likely an artifact due to this time step that I showed you before
Ole: in the
Ole: in the time like matrix.
So I,
Ole: although maybe there is also some consistent tracking I am. That would be so.
Ole: If I would be interested in that, I would. We could also try to
Ole: do the trf without
Ole: without 0 padding.
Ole: which means that instead of
Ole: instead of doing this doing this,
Ole: appending these zeros here, we will just cut off the stimulus at that point.
Ole: And this this way we can avoid this. the step function there. And let's see if we run it like that? Huh? Why did that not work?
Ole: Hmm.
that's weird.
Ole: Okay. That that should have worked. But I'm I'm not entirely sure I don't wanna
Ole: get into live debugging now
Ole: and then, if we we can also extend it more to the right, the time interval.
Ole: and we can see that towards the end. If we.
Ole: if you go like 2.8, or even like 1 s. make like a really long time interval from minus one to 1 s.
Ole: we can see that. Yeah, here, probably after like 400 ms or whatever, there's not too much too much change going on. So probably like 0 to.
Ole: or we can also probably something like 0 to
Ole: point 4 or whatever would be, would be a good choice here.
Ole: Any questions. So far.
Ole: you guys wanna do a little break? Or should we just power through the 2 h.
Stephanie Haro: one for break 0 for no break
Stephanie Haro: on the chat?
Stephanie Haro: No break is ringing
Ole: alright, but maybe then we should still take a a very quick, quick break to
Ole: account for the people who would like a little break, and II think I wouldn't mind either. So maybe
Ole: let's just take like the 4 min up to to and then meet back.
Stephanie Haro: Sounds good.
Ole: great.
Ole: all right. And let's
Ole: go on
Ole: I'll just skim over this section quickly, because I wanna continue with the rest of this stuff. So this is just to show you that there is a
Ole: the Trf class has like has a
Ole: method to export the data to the M. And E evoked class. And so this way, if you have M. And E installed you can. You can do that, and then you can access the
Ole: erez agmoni all the nice stuff for visualizations, and so on, that that M. And E. Offers. For example, what I'm doing here. I'm using this 2 M. And E evoked
Ole: class, which gives you out a list. and each element of the list is a evoked response. And like an er p
Ole: and you're getting one for each feature in the Trf. Here, in this case we only have a single feature, so we're only getting one. But it's still a list. So it's still a list with one element. So we have to take the 0 element of the list to get the actual evoked response. And then we can use this method plot joint here to get like a fancy plot that shows the the waveform and the scalp topographies at different points.
Ole: Okay. So now, to this question
Ole: of how to evaluate and optimize Tfs. And so we already played around a bit with the lambda parameter.
Ole: But
Ole: here we're gonna look at how to how to optimize it, how to how to choose the optimal lambda. But first let me show you
Ole: something else. Here I use the
Ole: I use the Trf. To predict method
Ole: with the Trf. We just trained to get a prediction of the stimulus and an estimate of the correlation between the prediction and the the response. So this is also so that if you call the Trf dot predict
Ole: method and you give it the stimulus and the response.
Ole: you'll get a prediction
Ole: and the correlation of the prediction, or with the the the actual response you can also call predict, using just the stimulus. But then you'll only get a prediction and no and no estimate of the accuracy.
Ole: and then we can compare it to the actual
Ole: to the actual stimulus at.
Ole: We can compare the response to the actual response and see that it has a correlation of 0 point 0 3 4. And this is the average correlation across channels.
Ole: so this is the
Ole: the observed eg. In blue, and then the predicted eg. In orange, and you can see it looks very different, but apparently it has some correlation.
Ole: but
Ole: that is this correlation is actually too high. So this is subject to overfitting, which is because we're training
Ole: and testing on the same data. So we trained on all of the stimulus. And we evaluated on all of the stimulus and response
Ole: Umhm. And so, because the
Ole: that fitting a T Rf. Finds the optimum mapping between the
Ole: between the stimulus and the response, it will also, it will
Ole: generally it will correspond to some of the brain processes that are going on. But it will also just fit some of the noise that's in the data. And so probably we're expecting that the if we apply that model to
Ole: some new unseen data, you probably get a lower accuracy than you. We would get for the training data. and to
Ole: to overcome this over this problem of overfitting.
Ole: we can use a cross validation. So this is a function in this stats module called Cross Val. And what cross validation does is it splits the data into test and training data
Ole: trains under training data tests, the model on the testing data
Ole: and then rotate us around so that each segment is used for testing once. And then so we and then we get
Ole: we get the average
Ole: correlation across all these all these segments. And what we're seeing here, it's like point 0 3 2, which is a bit lower than 0 point 0 3, 4
Ole: So what I would like to discuss now is.
Ole: why do you think it's we see, like a rather small change? Could you think of other circumstances. Where the difference between this cross validated correlation and the the raw correlation up here might be bigger like. Are there any any rules to this, what? What factors
Ole: affect the the amount of overfitting in a particular model?
Ole: Anyone has an idea like if we say we would like to increase the amount of overfitting
Ole: zoom.
Ole: So there are 2 main factors that that go into this
Ole: one is the amount of data that we have, and the other is the amount of parameters in our model. So
Ole: the more parameters, free parameters in our model, the more the better it is able to just fit random noise in the data, and on the other way around, the the less data we have, the more parameters per data point you have for fitting the data. So the amount of overfitting is kind of this function that depends on the the relationship between the amount of
Ole: parameters and the amount of data. Here we only have a single feature, just an envelope. And we have loads of data like, I think, like, it's like
Ole: almost 40 min or something of data. So
Ole: there is really not a lot of overfitting that can be done by the envelope here, just because there is only a single feature, and there is a and there's a lot of data, one thing which
Ole: we could try out as If we
Ole: if we include increase the number of time lags. I'm actually, actually, I haven't tried this. So we use like 0 to 1 s.
Ole: We're getting actually a higher correlation point O 3 9.
Ole: And when we do that non cross. Well, we get point 0 3 7. Okay, so this the number of time Lex doesn't really seem to affect the amount of overfitting very much.
Ole: So let's just get back to our
Ole: value of point 4.
Ole: Okay, so this is why why one has to use a cross validation. If you wanna get a correct estimate of your your accuracy. And th, this can be much more severe in in bigger models.
Ole: Hmm.
Ole: okay. And so now, to the question of optimizing lambda.
Ole: So before I talked about this problem of multi-col linearity, and how
Ole: regularization deals with that. But regularization also has another important purpose.
Ole: and namely, it prevents you to it prevents you from overfitting to your training data. And and this allows you to, and then to make better predictions about new unseen data.
Ole: So
Ole: so this means by regularizing, you're like, purposefully underfitting the the training data and with the with the aim of making better predictions of new data. So here I've simulated some data.
Ole: just a little twod example. So those are just 2 2 normally distributed variables that are correlated. So II have. I have the variable X. This is just 100 samples from the random normal distribution.
Ole: And then Y is also just 100 samples from the normal distribution, and plus point 6 times x, so the correlation of those variables should be should be
Ole: point 6
Ole: and then we can. We can sample random points, from the
Ole: from those data which I do here with this random choice function. So we pick 5 samples here.
Ole: and then we fit a. We fit a linear regression. Using the same methods that I just described. In the beginning.
Ole: we use the linear regression
Ole: to to fit
Ole: to fit those sample points and
Ole: and then see which regularization value actually gives us the closest approximation to the population trend here, I mean, in reality, usually we don't know the populate. We don't know the whole population. Right? We just we just know the data from our like few participants that we sample from the population. But here, in this simulation we know the actual population trend. So this is kind of the ground. Truth is this black dashed line, and then we can see which
Ole: which which
Ole: regularized regression comes closest to that. And here you can see that
Ole: with like
Ole: lambda equal one, maybe lambda equal. 2 or something would be even better. But this is the one that approximate, that gives you the best approximation
Ole: to the, to the trend based on our training data.
Ole: So now, what I would like to do, since all of this is random, we'll get a different result every time we rerun this cell.
Ole: So so I would like you to rerun it a few times in your own notebooks, and then look at under what circumstances you get. A low or high value of lambda will give you the best approximation of the of the population trend, and then you can also vary the number of samples you draw from the population.
Ole: So the this parameter here, and then see how that, how that affects the optimal value of lambda.
Ole: So I'll leave a quick break to
Ole: give you some time to try that out.
Ole: So do you think there is any rule that you can make out, whether, like in what kind of scenarios
Ole: a a low value of lambda gives you the best approximation of this black dashed line, and which scenarios a high value works better.
Stephanie Haro: 5 simulations seem to vary. Why, 3.
Stephanie Haro: From random sampling to random sampling.
Ole: Yeah.
Stephanie Haro: So it doesn't seem like I can learn anything
Stephanie Haro: a generalized rule.
Ole: So
Ole: if we try it out, for example, here, in this case, we're seeing that the that the
Ole: lambda equals 0 regression. So the unregularized regression seems to, or maybe like a very slight amount of regularization gives you the best
Ole: the best fit. And this is because here those 5 samples we picked, they give a very.
Ole: They're they're a very good representation of the overall population trend. And so, if that's the case, then of course, you don't need any regularization, because your your your sample, is very representative of all the other data. So based on this sample, I will be able to pretty. You get a pretty good prediction of all of the other data. Well, if I
Ole: rerun it, for example, like this. Now the now the trend, the the raw trend in my sample is is much steeper than in the actual
Ole: than in the actual population. And so in this case, I can. I can regularize, because what regularization does is it penalizes the the weights. So it kind of flattens out the line like, if you use
Ole: Max like infinitely large regularization, you'll just get a you'll just get a horizontal line.
Ole: and so I can flatten it out to get it closer to this population trend. And then, of course. it's in some cases
Ole: those actually all work pretty well.
Ole: Damn II can't get like a really bad one
Ole: sometimes. Of course, you're yeah. Maybe like this.
Ole: Yeah, like, sometimes you're, of course, doomed right? Because if you pick like
Ole: samples that are completely unrepresentative of the general population.
Ole: Then there is like no amount of statistics that can that can save you.
Ole: And this is why, it's important to to have sufficient sample size, both like on a subject level and population level. So on a subject level meaning recording sufficiently long, eg. From each subject.
Ole: and then, of course, recording sufficient number of subjects. So if we we can check that effect by by changing the number of
Ole: samples here, maybe from 5 to to 10, so double the number of samples.
Ole: And here generally, what we're what we'll see, or maybe let's do it even higher. Let's edit up to 20.
Ole: So now it should be that most of the time, because I'm sampling
Ole: a large number most of the time my samples should be should be representative of the the whole population.
Ole: and so most of the time I should get it. The the best estimate of the population trend with a with a small or regularization value or no regularization at all.
Ole: And this is something we're also seeing in Trfs, right? That is the that the more data you have the less regularization you need.
Ole: And just because
Ole: just because you're with more data, you're going to get, you're going to be able to
Ole: approximate the the general population trend better?
Ole: Does that make sense?
Ole: I'm assuming? Yes.
Ole: But also what we're seeing here, right is that you? You can't really a priori and decide
what?
Ole: What value of lambda works best. Rather, you kind of have to try it out. And this is exactly what we're also doing with trfs, I mean in in Trfs, we don't. I mean, in actual like data analyses. We don't know our like ground truth, population.
Ole: But what we can do is, we can kind of simulate that by just subsampling from our our sample, and then
Ole: and then fitting on one, and then
Ole: fitting the model on one part and testing it on the other, and doing that with multiple regularization values and see under which regularization we, we gonna get the best, the best out of sample prediction.
Ole: And this is what's done.
Ole: If you pass multiple regularization values to Trf train, which is what I'm doing here.
Ole: So if we run that Trf. Train will run the cross validation for each value in this in this array.
Ole: and then pick the best one to train the final model.
Ole: and then we can. We can plot the result. And we can see here in this case
Ole: the value of 26.8 3 was optimal, but the the values are spaced a bit weirdly because they're taking along a a log spaced
Ole: access.
Ole: But this is generally a good good strategy for for spacing the the values.
Ole: because the the higher the regularization is the the
Ole: the
Ole: the bigger the change in regularization has to be in order to have, like the same effect, so increasing from increasing from 100 to
Ole: for lambda, from one to 10, has a bigger effect than increasing it from 10 to 20, and so on.
Ole: And so now here, this means we're getting the best prediction with lambda equals 26.8 3. Again, the changes are pretty small. It's like the it only goes up from point O. 3 2 to point O, 3 7.
Ole: That is, again, because the envelope model is a very simple model that has few parameters, meaning it is very robust. And and we have a lot of data. So our results should not depend that strongly on regularization.
Ole: But still there is, there is definitely an effect.
Ole: And someone asked earlier, like, if, like, larger regularization is always better. And that is not the case. So here.
Ole: in this, in this example, and also here in the, in the, in this result, you can clearly see, there is like a certain optimum like in the beginning it gets better, and then it gets worse.
Ole: What you'll see most of the times is that the not not here? Because again, the envelope is such a simple model, but often you'll see like a steeper.
Ole: I mean, if we would extend the access, probably you would see it. Maybe, if we
Ole: do it like that.
Ole: You're usually seeing a more of a of a change towards the left side, so like in the sense that under regularizing the model has a a stronger effect than over regularizing it.
Ole: And I'm hoping. and I can. I can show that to you here.
Ole: No, apparently this model really doesn't need need much regularization, like often at at some point. If you go lower and lower at some point, you'll just see a complete drop off, because you'll end up with
Ole: with essentially this kind of model here. Right?
Ole: But I want to continue. Now, yeah, I have a question regarding the different subject. So you would try to optimize a lambda for each subject? Or would you try to optimize over all subjects? And how do you decide that?
Ole: So that depends on what you want to do.
Ole: so usually you would if you're like interested in the question of how strongly is that feature reflected in the brain response. So you're interested in that that accuracy estimate that that R value.
Ole: you would definitely do it on a subject by subject basis, because the data of some subjects
Ole: would be noisier than those of others, and then the noisier data might need more regularization.
Ole: But there is a problem that if you're interested in in directly looking at the weights. For example, I recently like
Ole: worked on a project where I was comparing the weights of the Trf. Across different features
Ole: and across subjects.
Ole: If you do that, it's problematic to use different regularization values. Because, as I showed you here. if you increase the regularization.
Ole: you're changing, you're drastically changing the scale of the y-axis and so
Ole: and and so it's hard to compare the the Trf. Weights of 2 different subjects if you use very different regularization. So in this case, I would actually try to find one value that works okay across all subjects and and keep it keep it the same.
Ole: And so it kind of depends on what you want to do in your analysis.
Camila Zugarramurdi: Sorry to follow up. How would you go about doing that manually? Would just pick up you and run it over different subjects and see what works best for all of them.
Ole: So so what I did in what I did in the in my project is, I just did the II optimized Lambda for each subject, and then just took the mean of the best
Ole: lambda. and then took that value for every subject.
Ole: So for for here we would get like a hundred now, because of like change, the the spacing of it. So apparently it arrives at a different value. And then for another subject, you get 1,000, and for another you get 10, and so on, and you just average them. And then you have, like the average best lambda. And you you using that for every subject.
Ole: And then it's probably not gonna be perfect. But it's it's gonna be okay for most subjects, at least.
Ole: is that, does it? Answer your question?
Camila Zugarramurdi: Yes, thanks.
Ole: And but actually, but this optimizing lambda introduces another source of overfitting.
Ole: because now, if we if we use that our value, we're getting out here
Ole: mit Ctl, and if we use that as our final estimate of the models. Accuracy that will be that will be inflated again. Because now, while the while we use cross-validation, meaning that the the Trf waits. 123
Ole: are only fitted on the train data. The lambda optimization uses both the train and the test data. So we're we have another free parameter
Ole: that we're fitting on the on the whole data set. So if we want to overcome that, we have to, we have to introduce a second level of of cross validation
Ole: which oh, wait! I I'm gonna like, reset my
Ole: range of regularization values.
Ole: and this can be done with this function called Nested Cross Val here, which is also the stats module.
Ole: and what this will do it will, it will do. 2 nested cross validation loops so in the in the outer loop. You split the data into a train and test set.
Ole: Then the train set goes into another cross validation loop where you fit a model, using each lambda, and then select the select the best lambda, and then that final model you'll test on the test set that you've
Ole: split out in the beginning.
Ole: And then you're gonna rotate through both loops, meaning for for each rotation of the outer cross validation loop. You're gonna do like a complete rotation of the inner cross validation loop.
Ole: And then and then you you're going to get one R value for each segment in the outer loop
Ole: and also one value for the best lambda.
Ole: And then the mean of this
Ole: of this list is A is an unbiased estimate of the model's accuracy, because the the test data
Ole: was never actually seen by the model fitting procedure. And is that is the idea of that clear? Why, why, there would be something you would want to do
Ole: okay? And but also here again, in this case, we see that it doesn't make a big difference again, because the I guess here, we're getting like something like point 3 7 as the best result. And and this is also what this function spits out. So. And this is again, because
Ole: we can see here the regularization doesn't have a really big effect right there. There. Most values we're testing are like, equally good. Almost so it it doesn't seem like we're doing a lot of overfitting here.
Ole: So and then to
Ole: finish that section, I wanna introduce the topic of statistical inference. maybe everyone start running this section because it takes a bit.
Ole: Now, we, we estimated this. Okay, here I had point 3 3. I think I used slightly different parameters here. It's like 0 point 0 3 7. But it doesn't really matter. So let's say, we've observed this
Ole: average correlation between the envelope and the response is
Ole: point O, 3, 3.
Ole: And now we wonder whether this is a high or a low value.
Ole: And so we want to do some statistical tests to assess whether whether this is a significant, this finding
Ole: and to do this, I already mentioned this in the talk we can use this function called permutation distribution to randomly permute the stimulus in the response and and get the correlation that happens just due to chance when there is no actual relationship between stimulus and response. I think this will run something like 3 min, or whatever.
Ole: So from this code we're getting, we're getting a histogram. or are we getting values that we can plot as a histogram where each where each
Ole: value represents one
Ole: random correlation, and then we can compare our actually observed value, which is here. The black dashed line to this distribution to assess the significance of our finding
Ole: So maybe, could you
Ole: take take a moment and think about
Ole: how you would now from this distribution. So here we have a hundred values for correlations. How can we calculate the p-value of our actual correlation from this from this list.
Ole: Anyone has an idea. I mean, it doesn't need to be. The full implementation can also just be like, conceptually.
Stephanie Haro: Camilla has a suggestion on the chat.
Ole: Yeah, yeah, this is exactly right. So the
Ole: so since the P value is just the probability of the
Ole: to observe the the data under the null hypothesis we can just count how often that value appears. And in python we can. We can just use
Ole: we can just use
Ole: logical indexing, so we can use our or a logical operation are permuted, and we can say where it is
Ole: greater or equal to wait. How did that? II think we said point O 3 3, or something like that.
Ole: And what that will give me is just a list
Ole: of Boolean values so either true or false, for every element.
Ole: and then I can. I can sum that.
and
Ole: And divide that by
Ole: by my number of permutations.
Ole: And this is the p-value, although there's a problem, because there's not a single observation, it's 0. But of course the p-value can't be 0. So in this case.
Ole: So we have to. We have to check that.
Ole: But if it's greater than 0 meaning we have some values.
Ole: or maybe we can, maybe we'll just.
Ole: We'll just call this N.
Ole: And if N is greater 0,
Ole: he is n
Ole: divided by the number of mutations, and if we and if we don't
Ole: we can make it a. We can make it a string.
Ole: and
Ole: else we can say that P is
Ole: smaller than
Ole: smaller.
Ole: and one
Ole: divided by N. Permute.
Ole: Why did this not work?
Ole: Oh, I forgot an equal sign here.
Ole: and that weekend
Ole: printout, P. So this means the number of permutations you're using, the
Ole: determines how specific your your test can be. So if I use the 1,000 permutations, I would be able to say, P, smaller point, O, one. But yeah, I guess from this plot. It's pretty obvious that the correlation is not random.
Ole: So now, okay, we still have about 25 min. Now, I wanna start to look into multivariate models.
Ole: So so far we've only used the envelope. Now, I want to introduce another feature
Ole: and to just like, keep it simple.
Ole: I'm here. I'm computing the acoustic envelope, the acoustic onsets which we get by computing the first derivative of the envelope
Ole: and then halfway rectifying that meaning. We're setting everything that is below 0 to 0. So we're only keeping the positive parts
Ole: erez agmoni. And if you, if we plot that we can see this is the onset, the orange line, and so we get
Ole: positive value is wherever there is a rising edge in the envelope. and this is known to be a good predictor of auditory brain activity are these, like sharp.
Ole: sharp increases in the acoustic input's power and the.
Ole: And then we can create a multivariate stimulus matrix by just
Ole: concatenating our, our exact existing stimuli.
Ole: and the rest works exactly the same. So you do trf dot train with your list of regularization values
Ole: and then plot the result.
Ole: I've had an exercise here, but I think for for time's sake we're going to skip that if that's okay for you, and maybe you can do it later if you want, and that would be to estimate the prediction, accuracy, using cross validation of this multivariate model, and then compare that to the, to the univariate envelope. Trf.
Ole: so potentially, this one should give us a better prediction because it's able to capture some a activity that the envelope model didn't. What you can see already is that the the the envelope? Weights are getting smaller because some of the variants that was formally explained by the the envelope, is now captured by this onset feature.
Ole: One thing that has to be considered with the with multivariate models is that you're you're only fitting one value of lambda for both features. And so it might. And so it might be that you end up with some sub optimal solution, because the the
Ole: the value for lambda that is optimal for one feature might not be optimal for another.
Ole: and one way to to get around that is to use what we call banded ridge regression.
Ole: and in banded ridge regression you split
Ole: the your features into different bands and then fit Lambda separately for each band. For example, if you had. if you had a 16 band spectrogram
Ole: and 39 phonemes, you could say the 16 band spectrogram is one band and the 39 phones are another band, and you're gonna optimize lambda separately for them. Here in our case both bands are just of size, one the envelope and the onsets, we consider a different band.
Ole: and then we can create our Trf instance and say, Method equals bandit.
Ole: And and this tells the Trf.
Ole: the Trf. That you want to use Banded Ridge regression.
Ole: And the the rest is exactly the same with the addition that you have to provide an argument called bands. And this is the this is the size of the feature bands in your stimulus that you wanna fit separate Lambdas to. So here in this case we're gonna say 1 one, because both of our feature bands have size, one
Ole: mit Ctl and the envelope and the onsets. If, for example, again, my my example from just now. If you had a 16 band spectrogram, and a 39 and 39 phonemes 2,
Ole: you would say 1639 here.
Ole: Also you'll you'll get an error message if the the sum of those is not equal to the size of your stimulus.
Ole: and if you run that this will actually run longer, because
Ole: we're doing an exhaustive search of all possible combination of lambda values, meaning that
Ole: meaning that the more different bands you're introducing the more different lambdas. You're fitting it. It. It increases the computation time exponentially. Because you're you're testing all possible combinations of lambdas.
Ole: So this should be used with with some caution.
Ole: Also, again, now you're introducing more free parameters. So also potentially more overfitting.
Ole: And then, if that has finished running, we can print out the result, and you'll see that there's actually a bit of a difference in the in the optimal lambda value estimated for both features.
Ole: Okay, I think I don't wanna necessarily wait on the result of this.
Ole: Are there any questions regarding Banded Ridge regression?
Ole: Alright, then I wanna come to this
Ole: to topic number 5, which is decoding.
Ole: So so far, we've predicted, eg. Responses from the envelope. But we can also turn this around and reconstruct the envelope from the brain responses. And that's th. That's what's called a backward model or a decoder
Ole: and this is
Ole: simply done by setting the direction parameter of the trf to minus one
Ole: and then just run the same model again with stimulus and response. And the trf will then just figure out that to use them the other way around the stimulus and the response.
Ole: this runs a lot longer because the model is a lot bigger.
Ole: Now, we're using now, our model has 128 features, because we have 128,
Ole: 128 channels. And if you, if you go
Ole: back all the way to the beginning
Ole: and look at the the auto covariance matrix.
Ole: So what you saw see is the auto covariance matrix has the has. The size is square, and the size is the number of time legs, times the number of features, because here's only a single feature. It's only the number, the number of time legs
Ole: erez agmoni. But now we have a 128 features. So we actually have 128
Ole: times 35 lags.
Ole: Oh, it's not. It's the
Ole: not running, because we're because it's still stuck on the other computation. But so this auto covariance matrix becomes a lot bigger. So our our matrix
Ole: computations are more demanding now. And this is this is why and the the backward model takes a lot longer. Also.
Ole: the free Google colab doesn't give us the best
Ole: the best CPU so this this takes for a bit.
Ole: But then, once it's once it's finished, we can also plot the result. And we can see now there is actually a much higher correlation
Ole: point one before it was like 0 point 0 3 something. You can. Still, if you plot the the reconstructed envelope versus the actual 100, there's still not too much, you can see just by. But the the correlation value is a lot higher.
Ole: However, okay, now it's finished.
Ole: however, this is also more subject to overfitting because of the higher number of parameters. So if you're gonna which I'm I'm not doing that here because it would take a lot of time if we gonna do the cross validation
Ole: on the the backward model. And O, also, I didn't optimize lambda now, but if we did the cross validation here, we would see a much bigger difference between the raw and the cross validated R than we saw with the envelope model. And that is because with all of those parameters, the the model can just do a much better job of of overfitting the noise in the in the training data.
Ole: Another
Ole: another important thing to consider when working with backwards models
Ole: is how to interpret and the the Trf weights.
Ole: So that is because
Ole: because the the backward Trf does not only amplify activity that is, of interest, but also suppresses activity that is not of interest. It kind of acts like a filter to extract this this latent stimulus factor from the data.
Ole: And this means that you can't. So in the forward model you can.
Ole: there is like a physiological interpretation. If we go back to maybe the Topopopo plots. There is like a physiological interpretation
Ole: of the weights, meaning that you can take this distribution of weights
Ole: to say, Well, if there's an increase
Ole: in the envelope
Ole: you're expecting at 130 ms. You're expecting positive deflections in the front of the scalp and negative deflections in the back of the scalp 100.
Ole: So so you can directly conclude from the model weight to like the the physiological activity. Basically, this is kind of like a
Ole: like a approximation of an evoked response potential.
Ole: But this interpretation is not valid when using a decoder, because the decoder does not have the aim of predicting predicting brain responses, but rather extracting the stimulus feature from the brain response. There is a
Ole: there's a really good, really good paper I linked here, that is, that deals is specifically about this issue on how to interpret or
Ole: why you can't interpret decoder weights just in the raw
Ole: but one thing you can do, which is a method introduced in that paper is, you can transform those backward weights to a forward model.
Ole: So there's a method that's called 2 forward. and that will give you the will, transform the backward models, weights
Ole: in a way that you can interpret them. In the same way you would interpret a forward model.
Ole: And we, when we compare that we see that the this, like transformed
Ole: envelope backward Trf is sort of a smooth out version of the the envelope. Trf, and this is due to these due to these filtering properties of the of the
Ole: the backward. T, rf.
Ole: so this is what you if you wanna interpret top topography plots of a of a backward T rf, this is usually what you should do is this this transformation to a forward model. And in this paper you can. You can read some some more detailed reasoning on on why? This is the case.
Ole: Any questions on that
Stephanie Haro: has that 2 forward functionality always been in the Matlab implementation of Mtrf. Toolbox, or there is a similar function there. I think there's something like Mtrf. Transform, or whatever.
Ole: So they they, they definitely also implement that method from this paper. I'm not. I'm not 100% sure how the function is called.
Ole: Also, it might be that it was added in the second in the second version of the Mtrf. Toolbox.
Stephanie Haro: Thanks.
Ole: Okay, cool. So I think now, I can actually, this is the last section, and we have, like, maybe 10 min left. So that works perfectly. Before I start talking, I I'll run this code because it also takes a while.
Ole: so for the last section, I want to address the question of how much data you need. So before we already saw that the amount of data determines
Ole: how good your correlations are, how much regularization you need, and so on. But usually, when you design an experiment.
you have to. You have to. A prior. You have to in the beginning decide on like, how much data you're gonna acquire from each subject.
Ole: So let's for now act as if we're in the process of piloting a new experiment.
Ole: And we wanna find out how how much data should we collect? And so one way to to answer this or address this question would be.
Ole: you just record a few subjects like maybe 2 or 3 people
Ole: record a lot of data from them. and then. like simulate
Ole: sort of new experiments from this data by splitting this data into short segments
Ole: and then rerunning your analysis with increasing amounts of data.
Ole: And then we can see how the how the prediction accuracy is is changing
Ole: as a function of the amount of data we're putting in and
Ole: This is what I'm doing here. So I decided on a segment duration of 5 s. Here you could also choose something else.
Ole: And then I'm
Ole: I'm splitting. I'm taking each the each trial, and I'm splitting it into 5 s long bits.
Ole: so that in the end I'll end up with with something like 700, something, five-second trials.
Ole: And then I can.
Ole: I can select varying amounts of training data.
Ole: Here I go up. This is in seconds. So I go up from 30 seconds of data up to the whole data. Set
Ole: and then for each.
Ole: for each duration, we'll do the whole thing like fit. A fit, a model and and get a prediction. Accuracy, or one thing is in right. At the very beginning I exclude the last 2 trials for testing. That is because we need, we just need a certain amount of of data to get a good estimate of our accuracy. Like, if I only have 30 s of data.
Ole: And and I just test with cross validation in the end, I'm I'm evaluating on only 5 s of data. And this would be.
Ole: And this is much too little to get a good to get a good estimate of the accuracy. And if you have only had, like very, very short testing data, you're gonna get very extreme results. So on on average, you you'll probably get you'll get the right result. But then in many, in many runs, you get like a very large negative number and other runs, you'll get a very large positive number, and you would have to rerun the analysis many, many times until you converge on the the actual mean value.
Ole: And so this is why I'm I'm excluding in the beginning 2 segments for testing, and then I'm testing all of those models on the same.
Ole: on the same segments.
Ole: and then what we can see
when we
Ole: when we plot the results. So here I plot, the correlation coefficient
Ole: for each training duration and the optimal value of lambda.
Ole: You can see that
Ole: the correlation is increasing. So here there's this first
Ole: step was a bit weird. But there is just some chance involved, because we're randomly sampling from the data. And so
Ole: if I would do such an analysis.
Ole: if I would actually do such an analysis I would rerun this whole pipeline like multiple times.
Ole: like maybe a hundred times, and then average. And then we wouldn't see, probably wouldn't see something like that. But now here you can nicely see that there seems to be like this sort of
Ole: logistic curve.
Ole: or like exponential curve here that that converges at some limit.
Ole: And here it seems like maybe after after
Ole: 1,000 or 1,500 s, you're you're not, or maybe even after 500 s. So if I would fit a curve, it would probably look something like this, 150.
Ole: So it seems like we're not. We're not getting
Ole: a lot of extra
Ole: power in our model. Just from
Ole: from including ever more data beyond this point. So here, I might actually conclude that okay getting
Ole: getting I don't know. 20 min of data from each subject is is fine, and that should be enough.
Ole: Of course this is only a single subject, and then what you can also see is that as the regular as the correlation is increasing, the value of lambda is is decreasing, and this, this is what we've seen before
Ole: that. With more data, you just need less regularization to to fit your model.
Ole: And do you have questions about that.
Ole: Should I look in the chat?
Ole: No, okay.
Ole: alright. And then if if some of you are interested in doing some more work with this data set and exploring it some more.
Ole: I have some bonus exercise here for you that you can do on your own
Ole: and that that is one thing
Ole: like, if you're in a situation where you can't acquire all the data you want. For example, you're recording from children, and you can't just record an hour of Eeg from children. This can be to some degree overcome by train, by by
Ole: training.
Ole: by recording a shorter amount of data from many subjects, and then training a Trf. For for each subject, and then averaging all the models, and then making predictions with these generalized models, and
Ole: this has been shown in different papers to to improve predictions, and so what one could do is
Ole: download the other subjects. So you would just use dataled again, and just, instead of the number of the subject put in an asterisk, and this will download all the subjects.
Ole: and then we can. We can go through all of the subjects in in the same
close, this
Ole: in the same fashion. We're doing it here. We can go through all of the subjects. And then for each training duration. Train a Trf. For the data of each subject.
and then
Ole: average all the Trfs.
Ole: the Trfs of all other subjects, meaning, if you want to predict the data of Subject 3. You're going to average the Trfs of all subjects. But subject 3,
Ole: and averaging is really simple with Mtrf Pi. So you can simply just
Ole: like you would with arrays. Just add add the Trfs and divide them. So an average T rf, would work like this. Trf, one plus 2 plus 3, and then divide it by 3,
Ole: and then predict the predict the data of that subject from the average Trf. Of all other subjects, and then you could compare the the the curve of accuracy over time between between average and specific models. And so what I would generally expect here is that the average model does better with limited amounts of training data. But then, if you have a lot of data, probably the
Ole: subject-specific model should outperform the the the generalized model, because it's just able to capture some, some subject, specific brain activity that is not present in the average.
Ole: So something for the interested. Who wanna explore this somewhat also? Yeah.
Camila Zugarramurdi: last question, sorry is this one of the cases where you would want to find the regularization parameter for all the subjects, because then you're gonna average them.
Yeah, yeah, exactly here. If you would, if you would have one subject that has a much higher lambda than the others. Then you would basically ignore that model and the average because the weights would be so low.
Mick Crosse: Yeah, very good point. I have a a question on that as well. Thanks very much for this. By the way, brilliant and great work on developing the toolbox.
Mick Crosse: So it's it's recommended not to actually average model weights when trying to combine subjects, or even trials, or whatever it is. And instead to actually sum over the covariance matrices?
Mick Crosse: And I'm wondering, is there a way using Mtf pi to do that? Could you import data?
Mick Crosse: For multiple subjects in order to generate one model.
Ole: II think there is not no so I mean, you can, of course, just use the the covariance functions from the toolbox. But
Ole: but that is actually a good idea to do that for the future.
Ole: Do you know a paper where they looked at the difference of averaging covariance matrices versus averaging
Ole: weights.
Mick Crosse: It was is a paper by Gustavus in 2,010
Mick Crosse: that showed that they're numerically different. there's a there's a kind of tutorial on
Mick Crosse: I think it's on auditory, specifically on Audrey. Attention decoding models that that also
Mick Crosse: talks a little bit about this. So I can throw a link into the chat if you want. Yeah, that would be great. I'll have a look at that.
Ole: I have to think a bit about how one would implement that to to use the covariances, because then you would need to pass in all of the data rather than just the the fitted. Trfs.
Ole: but I think it's potentially a useful feature. Yeah.
Mick Crosse: yeah, yeah, it is.
Mick Crosse: it. Also, it's kind of a nice way of
Mick Crosse: consolidating a lot of data as well without having the issues of you know, discontinuities, you know, between trials and things like that. It kind of circumvents all of that which is which is, could you also? With that, could you also have different regularization parameters per subject?
Ole: Or you're just well, that's the thing you. You wouldn't.
Mick Crosse: If you were to do it this way, you would have to
Mick Crosse: use an individual
Mick Crosse: or or just sorry, a single regularization problem which might be sub optimal in some ways.
Mick Crosse: Hmm.
Ole: alright.
Ole: Did you? Did you share the link to the paper?
Mick Crosse: I'll just put that in now
Ole: and then I think we should head back to the
Ole: the other meeting room.
Mick Crosse: Yeah.
Ole: Great. Thank you.
Mick Crosse: Okay.
Ole: Alright.
Stephanie Haro: Okay. Thank you so much for such a well formatted
Stephanie Haro: tutorial. There was a lot of embedded papers that people are going to find useful and lots of instructions. So I definitely will be going back to this file later on.
Stephanie Haro: and thanks for powering through Covid for us.
Ole: Yeah. Was a bit disappointed that the like second time Covid was almost as bad as the first time I thought it would be easier.
Stephanie Haro: Yeah.
Ole: well, is what is is.
Stephanie Haro: Well, thank you for following through even okay, so everyone
Stephanie Haro: if you wish head over to the other zoom link I could try to find it really quick and paste it on here for the QA.
Stephanie Haro: Thank you.
Ole: Alright, bye.