Scientist Stories: Merle Behr, Compositional Discrete Structures

Merle Behr presented methods for recovering hidden structures in mixed data. The first part covered finite-alphabet blind separation, used in cancer genetics and wireless communication, with an R package called haplosep. The second part introduced a random-forest approach to detect gene interactions, showing that true interactions have a specific prevalence signature that distinguishes them from noise.
Chapters
Yeah, good morning everyone. This is the Data Science Seminar. Welcome to Merle Baer, our speaker today. And now I will hand over to Jan Kleitzner to introduce Merle. It's your stage.
Thank you, and also a warm welcome from my side. It's my pleasure to introduce Merle Baer. She will be talking today, and I'm very much looking forward to the talk. Maybe to give a few sentences about her, she studied math, I believe, also in Göttingen, where she then received her PhD or worked on her PhD under the supervision of Alex Munk, or Munk, I guess, working on multi-scale change point methods and finite alphabet blind source separation. Today, she went, I think, to the United States with several scholarships and worked with Bin Yu also. And I saw a video of her talking online explaining her PhD thesis to laymen, so I'm confident that I might understand at least some parts of the talk today, but I'm still hoping that there is a lot of interesting math. And with that, I would like to hand over to Merle and hope Really looking very much forward to the talk.
Yeah, thank you very much for the kind introduction, and thank you very much for inviting me here today. So my presentation today will be about statistical recovery of compositional discrete structures, and this presentation will actually have two different parts which are kind of related, but still different. And the second part of this presentation is joint work with Professor Bin Yu and Chow Li and Yu Bang from UC Berkeley. So that is work which was initiated during my time as a postdoc at UC Berkeley. And the first part of the talk, that is about a follow-up project with my former PhD my supervisor, Professor Axel Munk from the University of Göttingen. So I want to start this presentation with sharing a quote with you, which I really like very much, and it's a quote from John Tukey, a very famous statistician, and he once said that the best thing about being a statistician is that you get to play in everyone's backyard. And for me, this captures really nicely what I really like about statistics, my field. So what are the different backyards that we will be looking at today? So maybe the first backyard statisticians like to call applications, but this is really what it is all about. We have some domain problem which we want to solve, and that's why we doing all our statistics and data analysis. And in the first part of the talk, so as I said, we'll have two different parts, and in the first part of the talk, the applications, we will look at three different applications or domain fields, and they are all quite different, but you will see that they share some underlying mathematical structure which we will then analyze. So the first application maybe fits best to this talk because it is about cancer, in particular cancer genetics. And then the second application is also from genetics but looks at data which comes from so-called evolve and resequence experiments. And the third application I want to quickly introduce you to is very different, namely from digital communication, but still they do share a common structure. And then in the second part of the talk, we will start by motivating the problem with an application again from genetics, namely higher-order gene interactions, so-called epistasis. So the second backyard, as a statistician and mathematician, we really want to go a step further and also understand the underlying mathematical structure which can describe this domain, the data from the domain field. And we want to analyze this theoretically. So the second backyard probably is theory.. And in the first part of the talk, we will look at so-called FABS structures, and we will look at identifiability and minimax optimal estimation procedures. And then the second part of the talk, we'll look at consistent Boolean interaction selection. And last but not least, the third, backyard computation. The best method is useless if we cannot compute it efficiently. And from the computational perspective, in the first part of the talk, an algorithm which will be important for us is a clustering algorithm, in particular a Lloyd's algorithm. And in the second part of the talk, what will be really crucial for us algorithmically is the random forest algorithm in order to compute things efficiently. So let's start, and as I said, we will first look at these three different applications, and you will notice a common structure for all of them. So the first one from cancer genetics particularly looks at so-called copy number variations. And as you all work for the German Cancer Research Center, at the German Cancer Research Center, you probably know this, or many of you know much more about this than I do. But what this is about in our context is that we have a tumor and, for example, this ugly purple guy here, and the tumor of a DNA is usually mutated. And one specific type of mutations are these so-called copy number variations. And that means that some parts of the DNA are either duplicated or deleted. So this is illustrated here in the middle. You see a healthy chromosome, and because we have one from our mother and one from our father, every region of the DNA will appear exactly 2 times. So for example, this red region here appears 2 times. However, if we have copy number variations, then some regions could be duplicated like here on the right. Now the red region appears 3 times or could also be deleted. So that are copy number variations. Unfortunately, things in practice do get more complicated. And the reason for that is if we take a probe of tumor cells, we usually do not just find one type of clones or type of tumor cells in there, but rather a few different types, so-called clones. And if we sequence this probe of cells with some bulk sequencing techniques, then what we actually observe is the mixture of these individual components. So what we're interested in are the individual clones, and of course there could also be some healthy cells among those in this probe. So we want to know something about the individual clones, but what we only get to see is a mixture of them. But there is one crucial information here, and this is by definition the copy number variations of the individual clones, because they correspond to actual numbers of copies of these pieces of DNA. They can only take values in a known finite alphabet, so they can only take the values 1, 2, 3, 4, 5, or maybe up to 8, but certainly it's a finite set. It's a relatively small set, and we'll see that this will be really important information for us to solve this problem. Problem. So we observe a mixture, but what we have is that we know that the individual components have this discrete finite alphabet structure. So let us look at the second example, a different data example. So the second example are so-called evolve and resequence experiments, and what people do there is they take some model organism, so for example this could be fruit flies,. And then they isolate these fruit flies, and then they apply some outer pressure on this population. And now this population is going to evolve. So these fruit flies, they will mate among each other. And because this population is locked up, no new fruit flies come in. This means that no new genetic information comes into this population. So what happens, these fruit flies fruit fly population becomes inbred over time, so only a few different types of haplotypes will survive in this population after a few generations. And now if we want to understand the genetic information in this population, it is really tedious to go ahead and sequence every single fruit fly in that population.. So what is usually done is instead we perform bulk sequencing. That means we take this whole group of fruit flies, we mix up their genetic information or material, and then we sequence that all together. Again, what we observe in this case is a mixture of these few different haplotypes which are present in this population. But what we are actually interested in are the individual haplotypes, so the individual types of genetic information which survived in this population. And again, we do have additional information, and this is that for the individual haplotypes, if we look at their SNP structure, then this can only take values in a given finite alphabet. So for a single haplotype, a SNP can either be 0 or 1, whereas in the bulk sequence we observe these mixed-up allele frequencies where all the different haplotypes contribute to. So again, we see a similar structure. On the one hand, what we observe is a mixture. On the other hand, the side information that we have is this finite alphabet of these individual components. So the third example, very different this time. So this time, the third example comes from wireless digital communication. So when we have our smartphones and we want to send a text message to a friend, then what happens is on our smartphone this message is encoded as a digital signal. So again, we see already digital, that is where the finite alphabet kicks in. But we now need to transfer this digital signal in a wireless way. So in order to do that, this digital signal needs to be transformed to an analog signal, and then it will be wirelessly transported through the air to some antenna which will receive this digital signal. So what happens on the antenna side is, for a single antenna, the different smartphones which all send these digital signals will all overlap, and all these different signals arrive there.. And we actually typically do not just have a single antenna but several antennas, and what we receive on the antenna side is then a mixture of these individual digital signals which came from the different smartphones. And we do observe several of those mixtures because we have several antennas at slightly different positions. So again, you can see that there's an analog structure here. On the one hand, what we observe or what the antenna observes is a mixture of different source components, and the additional side information that we have is that these source components came from a digital code which was binary. So the common theme of these three different examples that we introduced here is that we obtain a data matrix which corresponds to a linear mixture of finite alphabet sources with some mixing weights, and these mixing weights are continuous. So here we've written this down as a mathematical formula. So our observations, this is denoted as Y_ij. It has two indices because it is a matrix, a matrix with small n rows and capital M columns, and the individual components in this data matrix correspond to such a mixture where these source components are these finite alphabet components. So for example, for binary alphabet, this would mean that these elements here can only take the values 0 and 1. And they are mixed up with some mixing weights. And here in this presentation, we will assume that those mixing weights are positive and they sum up to 1. So this is, for example, the case in the cancer genetics example where those mixing weights correspond to the physical mixing proportions of the different clones in the tumor. It is also the case in the evolve and resequence experiments where again those mixing weights correspond to the fraction of the respective haplotypes haplotype in the generation under consideration. And we also need to assume that we do not observe this mixture exactly, but that there is some additional noise component to it. So let's look again at our three examples. What does this exactly correspond to? So in the cancer example, the small m, so that is what the sum here goes over, the different source components, that would correspond to the different tumor clones. So we would have small m different clones in our probe. The small n, that would correspond to the number of genome locations where either we have data for our copy number variations or also, yeah, or that would be the case in our tumor example. And the capital M, the number of mixtures, that in the tumor example would, for example, correspond to the number of different tissue samples that was obtained. So if we obtain tissue samples at different locations in the body or at different time points where we would expect a different fraction for the respective tumor clones. Similar in the above and resequence experiments, the small m, so the number of sources, that will correspond to the number of major haplotypes which dominate this population. The small n that corresponds again to the number of genome locations, so the SNPs where we get measured allele frequencies at, and the capital M, that would correspond to the number of generations. So remember, this population evolves over time and we would have different generations, we would have different proportions of these haplotypes in the generation. In the wireless digital communication example, the small m would correspond to the number of senders. So for example, that could be the number of mobile devices. The small n would correspond to the length of the code, and the capital M would correspond to the number of receiver antennas. So this is now this common data structure that we have at hand, and what we obtain is our data matrix. And what we are interested in are these finite alphabet sources as well as these mixing weights. So that is what we would like to recover from our data matrix. So as a first step to understand this general structure of this problem, it is really illustrative to start with a simple example. And as a simple example, let us look at the situation where we would have only 3 source functions. This could correspond in these evolve and resequence experiments to 3 different haplotypes. And we would have 2 mixtures, so that could correspond to 2 different generations at which we observe different fractions of these haplotypes in our population. And for the evolve and resequence experiments, the alphabet would be binary, so it would be the values 0 or 1. And then as another simplification to understand the structure, let us also assume that there would be no noise at all. So we would have no noise and we would directly observe this mixture of these finite alphabet sources, and that would be our data matrix. So these green dots, that could be now our observations, so the data that we observed at different genome locations. So that would be one genome location, another dot would be another genome location.. And the x and the y-axis here, that would correspond to the two different generations. So remember, because we have two different sources, an individual row in our data matrix has two dimensions, and these are the two dimensions that you see here on the x and the y-axis, and that could correspond to two different generations. And then this could be what we observe, these green dots.. So now what we want to know is how can we obtain from these observations these finite alphabet sources and the mixing ones. So the first observation that we can make is that the smallest possible value of observations that we could ever see in terms of its absolute norm would appear when all the source functions take the smallest alphabet value, so when they are all 0. And in that case, what we observe is the all-zero vector. So that's exactly the dot here in the bottom left, and we know that it must correspond to the situation where all the different haplotypes have a SNP value equal to 0, so our 0, 0, 0 vector. Okay, so we have observed the first, or recovered the first source component. Now we can go on and ask, well, what is about the observation with the second smallest absolute norm? And we realize that this has to be the one where again all the source functions are zero, but only the one source function with the smallest mixing components, again in terms of absolute norm, is one. And at that position, we will observe the source vector 1, 0, 0, and this observation exactly corresponds to the first mixing component w1. So now we have the first mixing component and we can kind of subtract that from our data points and look at what among the remaining observations is the smallest possible or the smallest observation in terms of absolute norm. And so that means we continue recursively. We again, we look at the smallest one, we realize okay, that must be w2, so we're now only the source function with the second smallest mixing weight is 1, and so on. And then in a recursive manner, we can work us through these individual data points and always make up what has to be the respective combination of source functions. And it turns out that the same thing can be done for an arbitrary alphabet, an arbitrary number of source functions, and an arbitrary number of mixtures. So we can always go ahead and make use of this discrete structure, and by doing that, recursively can work us through the individual data points by ordering with respect to their absolute norm and use that to reconstruct the individual combinations of discrete alphabet source function values. So what can go wrong with this? So I said that this always works, but it's actually not true. There are situations where this does not work, and it's important to figure out what exactly those situations are because that will tell us what we need in order to be identifiable. So if we think about this in more detail, what we realize is the first situation where, or one situation where we get a problem with this recursive algorithm is actually when two of these observations are really, really close to each other. So for example, if these two red observations are very close to each other such that they are basically the same. So when they are the same, then of course we have no chance to distinguish whether at that particular observation we observe the value 1, 0, 0 for the finite alphabet component or the value 0, 1, 0 for the finite alphabet component. So a first necessary condition in order to guarantee that this whole thing is actually identifiable is that different alphabet values indeed result in different mixture values. Because if that's not the case, then we have no chance to distinguish between these finite alphabet source values. So we can write that down in mathematical terms, and we just write down the minimum distance between these different mixture values, and that has to be larger than zero. There's a second situation where this recursive algorithm does not work, and that is if we don't observe those weight vectors at all on their own. So recall that we were really relying on going through this recursively and recovering one observation by the other one. However, if some of these different combinations are just not present in our data, then we cannot do that. And it turns out that a sufficient condition now for identifiability is to assume that these individual weight vectors, so the ones where exactly one of those source functions is active and the other ones are zero, is present among our data. That indeed is a sufficient condition for identifiability. Okay, so we have now worked out two different identifiability conditions, and if you want to look at precisely the proof why indeed they are sufficient for identifiability for an arbitrary alphabet, arbitrary number of source functions, arbitrary number of mixtures, you can look that up here in this paper. But now we really have the basis to go one step further and to become more practical, because recall that so far so far we were assuming that there's no noise at all, and this of course is not a realistic assumption. In practice, we will need to deal with a situation where our data is noisy. So in that case, what one can do is one can follow a very simple 2-step estimation procedure. So the first step of this estimation procedure is just a very simple clustering step where you take your observations at different genome locations and you apply a Lloyd's clustering to obtain mixture centres. And then the second step is you take those mixture centres and then you apply this recursive algorithm that we just explained, and then you use this to recover the source functions as well as the mixing weights. So what we can now show is if we take this very simple two-step algorithm and we also make a few additional assumptions, so for example, we assume that the noise is IID Gaussians, and we do assume basically a slight variation of these identifiability conditions that I've shown you in the previous slide, then we can show asymptotically as the sample sizes increases, so the number of observations as well as the number of mixtures go to infinity at a sufficient speed, then with high probability On the one hand, it will hold that we do indeed recover these finite alphabet sources exactly with this algorithm. And on the other hand, we do get an estimate for these mixing weights which has a 1 over n estimation rate. So if you take the square root on both sides, then you see this as a classical 1 over square root of n estimation rate. And this is quite interesting because what this tells us is, recall that that if actually those source functions would be known from the beginning, then this whole problem of this finite alphabet blind separation, that would just be linear regression. And in linear regression, we would know that under this Gaussian noise assumption, indeed, we would have the standard 1 over square root of n estimation rate. And what we see here is that even though those finite alphabet sources are not known, we still get this same parametric 1 over square root of estimation rate here. So we don't lose anything by not knowing this finite alphabet source function structure. We also have more applied work on this, on these finite alphabet blind separations, in particular in this context of Evolve and resequence experiments. There's also an R package which is called haplosep, and it's essentially an extension of this very simple algorithm, which particularly is tailored for the haplotype separation. You can read more about this in this recent publication here, and we found it really to work very well. It's a very fast algorithm. As you can see, there's really only two very simple components, and it does provide pretty accurate estimates for these different haplotypes and also for their weights, which means that you can actually see how these different dominant haplotypes evolve over time as the generations evolve. Okay, so now we come to the second part of the talk, and at first glance you will think it looks pretty different, but it's actually quite related to the first part of the talk. So in the second part of the talk, we will look at recovery of feature interactions with random forest. And we start off again by looking at a particular domain problem, particular application, namely the application of epistasis and gene interactions. So originally epistasis meant that we have some genes where one gene is masking the effect of the other gene. So the classical textbook example is where we look at the fur color of mice. We would have one gene which is responsible for pigmentation. So one allele, which would be the dominant allele, implies that the mouse can have pigmentation, whereas the other allele, which would be the recessive allele, implies that the mouse cannot have any pigmentation at all. And then we also have another gene which is responsible for the color of the fur. And there we would have, say, one allele which corresponds to a dark color, whereas the other allele corresponds to a light color. And in this situation, what we see is— so you see a table here where the rows correspond to the alleles which are inherited by the mother, and the columns correspond to those inherited by the father. And the different entries in this matrix correspond to what the offspring mouse then gets. And what we see here is that actually the first gene, the gene A, can mask the effect of the second gene, because if it inherits two of the recessive alleles, then the mouse has no pigmentation at all. And no matter what the other gene then looks like, we we cannot see that from the mouse. So the mouse just has no pigmentation, no matter what the other gene is, whether it corresponds to dark brown or light brown, we just don't see that in the phenotype. So that is a very classical textbook example of epistasis. But typically, genetics is much more complicated than this very simple relationship between phenotype and genotype. And this is why later on people then started to define epistasis more generally, particularly going back to a definition by Fisher. It is defined as the departure from additivity of genetic effects. So we again have some phenotype of interest. So that could be the fur color of a mouse, it could be our own hair color or whatever we're interested in. And what we want to find is genes which do have some nonlinear, non-additive effect with respect to this phenotype.. So that's now our task or our domain questions: how can we define or find these gene interactions? And what we will do here in the following, we will consider the random forest algorithm to help us with that. And therefore I'd like to quickly have a quick recap of the random forest algorithm just to make sure we're all on the same page.. So the random forest algorithm as input data is a classical supervised learning setting where we have labeled data. So we have pairs of data where we have a feature vector which here will be of length small p and we have a label which will be an element in R. And the random forest algorithm, what it does in order to fit this relationship between feature vector and label is it builds on decision trees. So the individual components of the random forest are decision trees, and these decision trees are grown using the CART algorithm. So the CART algorithm proceeds as follows: we start at the root node, so we proceed recursively, and as a first step, what we do at every node is we randomly select a subset of features. So we will not consider all p features, but only a subset of them. And then, given this randomly selected subset, we will select among those the feature which provides the best split of the data. So the feature— the data is then split into two halves. For example, we select feature x1 and we split it at a threshold 0.2, so So all the data which satisfies this inequality will go on the left side of the tree and the remaining data will go on the right side of the tree. So that introduces some hyperplane which provides a split into two groups of the data. And then we continue and at some point we have to stop and in the random forest algorithm we will grow those trees to purity. So we will only stop if we really cannot split any further.. So that is, for instance, when we have no observations left or just a single observation which cannot be split into two parts, or if all the observations which are left are already completely pure. So here, as a purity measure, we'll consider the variance. So that would mean that they all have zero variance. So then there's also no way to split them any further by gaining any decrease in impurity. And then the final prediction in a card is just the average of the responses which fall into that respective leaf node. Now the random forest algorithm does a very simple thing. It will just grow various of those decision trees, and then the overall prediction of the random forest algorithm is the average of the predictions of the cards.. So the classical implementation of random forest will also assume that every decision tree just sees a bootstrap subsample of the data. We will not consider that here, but if we have many data, it usually also doesn't make a too big difference. But still, what is important here for you to remember is that these decision trees will all be different, and the reason for that is because we had this randomness in selecting this random subset of features. So at every node, we only We're selecting a random subset of features, and this is why these decision trees will all be different and random indeed. So what does this have to do with interaction? So why do we consider the random forest when we want to extract interactions between features? So the reason why we do that is because the random forest algorithm, or the decision tree in particular, comes with a very natural interpretation of interaction. And this is why it's individual decision path in the tree. So if we look at such a tree and then we look at a particular path in the tree, then what this corresponds to is the following. So when feature x1 is smaller than something and feature x2 is larger than something, then we have this particular phenotype. So you see that there is just this very natural interpretation that features which appear together on such a decision path do correspond to an interaction. And this is exactly the key idea that we want to make use of in the random forest in order to extract interactions. And there has been various literature before which was exactly making use of this particular key idea. So this idea that features which often appear together in decision paths which were fit to this labeled data, they do correspond to interactions. And it was found that this seemed to perform pretty well in practice. So also for epistasis, this was applied previously, but also in particular in other genetics contexts, and it was found to provide reasonable empirical interactions. We also have a recent work where we suggests a particular pipeline for this problem of detecting epistasis, which was a collaboration with the John Zuckerberg Biohub initiative, and you can read more about this in this recent preprint. There's also software for this, so that's our IP3 pipeline. And one highlight about this which I would like to mention is that it takes SNP data as input, but it also provides you with and not just interactions on the SNP level, but also on the gene level. So this is evidence, empirical evidence, so to say, that this seemed to work in practice. But what was really missing so far was any theoretical evidence. So we can just do this, and it just does seem to work because we can validate it in some settings where we know some underlying biology about interactions, and it does just seem to recover the right thing. But what we didn't really understand so far was why does that really work and can we actually prove that it works? And this is something I'd like to give you some insights on now on the following slides. And as always, if we want to prove that something works, first we need to define under what conditions on the data do we want to show that it works. So we need to introduce some model, and the model we'll be working on is the locally spiky sparse model, which we introduce here. And it has the advantage that compared to other interaction models, it does not rely on any Lipschitz assumptions, so it does not rely on any continuity assumptions. Instead, it really considers a non-differentiable continuous Boolean interaction model where we have j different interactions. So here is our response for our labels, and given our features, we assume that it corresponds to a linear combination of j different interactions, and these individual interaction terms are additive Boolean functions.. So under this model, what we would like to recover are exactly these sets of features which interact with each other in this model. So you see a very simple example here. So say we would only have a single interaction, which would be between the first and the second feature, then how this could look like would be as follows. So you see in the response surface, we have this area which corresponds to some hyperplane with respect to these two different features. Where we do get an additive component towards this response function. We do also make a few further assumptions which I want to mention here. So as usual, we assume our data to be IID from some underlying distributions. We do assume that the beta and the gamma coefficients here in this model are bounded, and we assume the features to be uniformly distributed on this hypercube. Remember that for random forest in general, it is invariant under under monotone transformation, so you could generalize that just to features being independent with some density. We assume sub-Gaussian noise, and what we also assume is that those interaction sets are disjoint. So this is actually something which we need to assume for identifiability reasons, because otherwise this type of model just could happen not to be identifiable at all. And here I also quickly want to mention that although it looks different at first glance, it actually turns out that this particular type of model is very related to this FAPS model that we've seen previously. So indeed, you can rewrite this in such a form that it does exactly correspond to a FAPS model. However, in the FAPS setting, we were looking at an unsupervised setting. We didn't have labels, whereas here we are in a supervised setting where we have labelled data, so we have much more structure in this LSS model than we had in the FABS model. So before we, or when we now want to quantify this idea of features which often appear together on the decision paths interact with each other, we really need a measure of what it means that features often appear together. So we somehow need to mathematically quantify what do we mean by this phrase, features often appear together on the decision path. What does often mean? And in order to understand this or to do this, it really is important to understand what are the different sources of randomness that are present in the random forest algorithm. So the first source of randomness that we have is the data randomness. So we do have our data and we have one realization of this data, and this randomness is kind of unavoidable. We just have the data that we have and we cannot go in a parallel universe and collect more and more realizations of this data. So that randomness is unavoidable, but we do have additional sources of randomness in the random forest, and the first source of randomness is the randomness of the trees. So remember that in these decision trees, at every node we were randomly selecting a subset of features, and that source of randomness is something which we have at hand. And we also introduced a third source of randomness, and this is selecting randomly a path in this decision tree. So we do this in the following way: we start at the root of the decision tree and then we randomly either go left or right, and that will provide a randomness which creates a path that is actually weighted with the inverse of the length of the path. Okay, so, and these two sources of randomness are randomness introduced by the algorithm, and that means that we can generate as many realizations of this as we like., because we can just generate as many trees in the random forest as we like, so we can average out over these two types of randomness. And this is exactly what we do and how we quantify this features often appear together on the decision path. So we define this so-called depth-weighted prevalence, and for a given set of features, this is exactly the probability and probability now means we are averaging over the two sources of randomness which we generated ourselves, that these features appear together on such a randomly selected path in the decision tree conditioned on the data. Okay, so now we have this thing which we can just compute given the tree ensemble from a random forest, we can just compute it. And now comes the theorem that we've shown.. So if we take this depth-weighted prevalence, which we computed from the random forest tree ensemble, then we could show that under the LSS model, if we indeed have a true interaction, then this depth-weighted prevalence equals 2 to the power minus the power of the size of the interaction. However, if we do not have an interaction, if it's not an interaction, then we do have this upper bound where we see that the steps-weighted prevalence is indeed strictly smaller than that by some constant c, for which we also have an explicit expression. And based on the theorem, we now do see that we indeed get consistent recovery of these Boolean interaction components, because all we need to do is we generate the tree ensemble from the random forest, and then we just select those sets of features where the steps-weighted prevalence does best guess correspond to 2 to the minus power of the size of the interaction. And that gives us consistent recovery of these Boolean interactions. So we indeed have shown that what seems to work empirically indeed has a theoretical basis with provable consistent recovery of these Boolean interactions. So I'm almost at the end. I just want to mention two additional insights that we gain from this analysis, from the theoretical analysis. So first of all, maybe the most obvious insight that we see is this term here decreases exponentially with the size of the interaction. So that really shows us when we want to recover higher-order interaction, it becomes exponentially more difficult. And the way you can think of it is is that these Boolean interactions, the size where they are active, that size of this hypercube will shrink down exponentially as we have more features involved in that interaction. So in order to actually see data which fall into that hypercube, we will need to sample exponentially more data so that we get a chance that some observation falls into that area where this interaction is active. So that's one thing, one takeaway. So we really can quantify with this that this scales exponentially with the size of the interaction. And then there's another insight I want to highlight, which is more practical, and that is that our analysis also provides a suggestion for an optimal sampling fraction in the random forest. So this sampling fraction is this parameter p. Usually call m-try. And it turns out that this constant c here in this theorem actually depends on that m-try parameter, and it also depends on the overall sparsity of the model, so the overall number of signal features in the model. And we can work out what the optimal constant is, and the optimal constant of course would be the one where this is as big as possible, because then our gap between true interactions and false interactions becomes as big as possible. And that happens when the m try parameter is of this form. And with this, we can now recover some of the default implementations for random forest, namely if the number of signal features equals a third of the overall features. So if that means two-thirds of our features are noisy features and only one-third of the features are signal features, then what we get for m try is this p over 3, which is one of the default parameters in many implementations. However, when we're in a situation where we are very sparse, so where we have much, much more noisy features than signal features, then we see that the m_tri parameter which would be optimal with respect to this interaction selection guarantee would actually be larger and it would be equal to p/2. So that maybe gives us a hint that if we're interested in interaction selection, we would rather want to select an m_tri parameter which is larger than the default setting in the situation where we believe that our model is very sparse. So I'm at the end of my presentation. I'd like to summarize. So in the first part of the talk, we started with multivariate finite alphabet blind separation. And among others, we motivated this with these Evolve and ReSequence experiments. And here you can find, again, the ref where we have a detailed applied analysis for this particular type of algorithm in these E&R contexts or for haplotype separation more generally. And here is the reference for the theoretical analysis I talked about. It's a preprint, and I also want to mention that you can also find a data example from cancer genetics in the this preprint here. I also want to mention two other related references here. So one which particularly is about identifiability in this model, and the other one is more on this change point estimation problem, which really where typical applications would be these copy number variations where you have this change point structure as well.. And then the second part of the talk, we started by motivating the structure we looked at with detecting epistasis and higher or higher order gene interactions, where we have this epitree pipeline, which you can find in this preprint here. And then we looked— we really wanted to complement this empirical evidence with some theoretical evidence that indeed shows that this type of reasoning does provide consistent interaction selection, and that you can find here in this recent preprint.
And yeah, with this, thank you very much for your attention. Thank you very much for this excellent talk and the nice explanations. I hope people have some interesting questions for you as well. So please just ask them in the chat or raise your hand if you have any. And maybe wait one second. There's not— none yet. So maybe I'll start with one that is related to the first part and the identifiability of the FABS problem. I wonder how does the identifiability depend on the mixing coefficients?
I assume for simulations you would probably use quite similar-ish or—
Yeah, so let me go to that slide.
I think it's here. You mentioned when they're quite equal, then you have a problem? Exactly. So this, the mixing, so we had two conditions here and really these two conditions, the first one is a condition on the mixing weights, whereas the second one is a condition on the source functions. So how does this depend on the mixing weights? So that's when we have to look at the first condition.. So exactly when they are very similar, this is a property of the mixing weights. So for example, if two mixing weights would be exactly the same, that would be a situation where these two different weights would be the same. And of course, in practice, you can be very unlucky and just run into a situation— I mean, you can numerically write down mixing weights ways where this is zero. So in principle, you could be very unlucky that it's just exactly zero, and then you have no chance to identify these individual components because they would just not be identifiable. However, it turns out that if you assume that your weights are in a sense random, so for example, that they would be uniformly distributed between 0 and 1, then as the number of mixtures, so this capital M as increases very, very quickly, it becomes extremely unlikely that you do have a situation where this is zero. So actually, the probability, if you have uniform weights, the probability that this is zero is exactly zero. And in particular, if the number of mixing weights increases, so if you have, for example, more generations or if you have more samples from your tumor, then the probability that this happens, even that this is larger than some threshold, this is exactly what you will then need precisely in the theorem. So 0 is not enough, you need it to be bounded away from 0 actually.
But that will also hold true with a high probability as you get more and more mixtures.
And these conditions are now just for the case without noise? Right, because they are conditions just on the mixing weights and on the source function. So they are just conditions on the underlying structure of the problem, and then on top of that, if we actually want to recover the this from data, then of course we also need to take into account that there is noise. And the conditions on the noise, and for this particular theorem that was then on the next slide, was that the noise is Gaussian with some variance component which does not have to be known, but which should be the same for the individual noise dimensions.
This depends on the constant that is There are zeros. Say it again, the question?
This maximum variance for the noise depends on the constant that you have in the— Oh, right, yeah, so you're then wondering where do these constants enter. Exactly, so the variance in the noise here, that will actually enter in a constant which is here hidden by this twiggles and it also enters in a constant which is hidden in this "with high probability" because "with high probability" actually has an exact expression, and in that expression, there will also appear the noise variance. So you're absolutely right.
As this noise variance increases, this probability will also decrease, and it influences the constants which appear in the theorem. Yeah, thank you. I think you explained everything very well. Since there so far no questions, I'd again encourage you to ask any questions or raise hands. And in the meantime, I have one more, if you, if you have. So does your algorithm need to know the number of components?
So for example, the number of clones in the Yeah, that's a good point. So actually, yeah, so this one here, as it's written down here, does require to know this. A practical way to estimate this in practice is to look at a singular value decomposition of your data matrix. And the reason for that is under these identifiability conditions, which were basically saying that we have enough variations in our source function you can— let's see if I— you can show that the noiseless data matrix— where do I have it? Yeah, so that this matrix here, that the rank of that matrix is equal to the number of source functions. So now what you can do is you could just look at the problem of having a low rank, a noisy low rank matrix, and estimating the rank of the rank of that noisy low-rank matrix, and that, for example, you can do by looking at the singular values of your data matrix. And in our R package for HabloSep, there's also an implementation of that available, but we did not show any theory about that, so in our theorem here, we would assume the correct number of source components to be provided, which can be estimated via the singular values.
Yeah. Okay, thank you. There are no further questions. I'd like to thank you again and ask everyone to thank you as well and hope you might reconsider publishing research again.
I think you switched to FAMR. So, well, thanks for the talk and looking forward. Yeah, thank you very much for having me. Yeah, thank you also from my side.
So thank you for being here.
Thank you also to the audience.
Goodbye.
Bye, thank you.
More from Joshua

Scientist Stories: Hans Waterham, The value of phenotypic assays in peroxisomal disorders

Scientist Stories: Joel Zlotogora, Genetic carrier screening programs for reproductive purposes








