Custom bootstrapping scheme

[tecosaur]

Hello,

I'm interested in applying EvoTrees to some of my work, however there's a particular constraint that I need to apply to my dataset that I can't see any easy way to incorporate into EvoTrees.

I'm working with a dataset with many "clumps" of rows, which are highly correlated with each other. When fitting a plain random forest to this data previously we adjusted for this by employing a stratified sampling technique: each round of training we randomly selected one row from each "clump". This preserves the assumption of independence during bootstrapping.

If I'm mistaken and there is a way to acomplish this with EvoTrees (even if it means reaching into internals a bit) please let me know how I might do this, otherwise this is a feature (a custom/stratified bootstrap scheme) that would be tremendously useful to have (for my use case at least :).

[jeremiedb]

Thanks for your inquiry.

I tihnk it would be effectively tricky to emulate this sampling behavior directly thorugh current EvoTrees' API.

That being said, an adhoc hacking may not be too cumbersome. The sampling operation, performed at each iteration round, is done here:

EvoTrees.jl/src/fit.jl

Line 11 in 2435c1c

cache.nodes[1].is = subsample(cache.left, cache.is, cache.mask_cond, params.rowsample, params.rng)

The implementation of subsample isn't that trivial, but only due to performance optimisation reason.
The only requirement for subsample is the return a view of observations ids (a vector of 1:nobs).

EvoTrees.jl/src/subsample.jl

Line 6 in 2435c1c

function subsample(left::AbstractVector, is::AbstractVector, mask_cond::AbstractVector{UInt8}, rowsample::AbstractFloat, rng)

The challenge here would be to make available to a custom subsample the information about the clump/group each id belongs too. This could be done eventually cleanly by add new information variable kwargs, similar to current weight_name/offset_varname (ex: group_name), then add this grouping vector to the cache at model initialisation step.

I'd note that I see 2 different expectation about the handling of groups in tabular data:

1. subsample N observations from each group (what you're referring to here, where N=1)
2. subsample full set of observations belonging to N groups (ex: in ranking task, such subsampling is useful if one want to build a loss function over all the observation associated to a query).

For your use case, I'm not clear if the non-dependency of your sampling would actually be an issue. Have you compared the out-of-sample performance if using vanilla subsampling? If your groups have largely varying size and the concern is that some groups weight may be over/under represented, an approach would be to use the existing weight_name to use it to specify weights inversely proportional to the clump over-representation. (ex: group A has 100 obs and B 300 obs, then weights could be set respectively set to 0.75 and 0.25 to each member of A and B.

[tecosaur]

Thanks for the detailed response Jeremie!

To help ground the discussion in a specific example, and avoid the XY problem, the specific case I'm working with is where I have a dataset of gene variants (with one or more variants per gene in the dataset), but most of the feature columns are per-gene, not per-variant. I'm wanting to build models that learn to distinguish between two classes of variants, one of which is much larger (up to 10,000x in some cases) than the other.

I haven't directly compared vanilla/weighted sub-sampling. While the weighting is a factor, the key concern here is that variants of the same gene (which share feature information) will be split across the test/train sets, violating the IID rows assumption of bootstrapping/random forest construction. Without IID rows, the OOB evaluation will be overly optimistic. This also introduces (in the case of random forests) or increases (with boosted RF) the risk of over-fitting and complicates downstream statistical analysis.

To account for these factors, I've ended up with this custom bootstrap scheme:

1. (prep) I split the data into two halves: one for each category
2. (prep) I split each category into genes, as a vector of UnitRange{Int} row indices into the data matrix
3. (train set) I perform a sample of genes from the smaller category, and sample from min(50 * length(smaller), length(larger)) genes of the larger category
4. (train set) Select one variant at random from each gene
5. (test set) Use all variants from genes not selected to be in the training set

From the code that you pointed to, it looks like it might be possible to adapt subsample and grow_evotree! to take an extra keyword arguments or two to provide the vector of UnitRange{Int} row indices I mentioned and the category split.

How does this all sound to you?

[jeremiedb]

Is the following a proper depiction of the nature of the dataset:

 Row │ geneID  variantID  x1         y
     │ String  Int64      Float64    Bool
─────┼─────────────────────────────────────
   1 │ A               1  0.675536   false
   2 │ A               2  0.0461865   true
   3 │ A               3  0.0104879   true
   4 │ B               1  0.57438    false
   5 │ B               2  0.523138   false

Is there some inherent correlation between the number of of variants and the prevalence of the variant class? Ie, is
E[Y | number_of_variants] = E[Y]?

If the above condition is respected, my intuition would be that simply splitting the train-eval based on geneID would be sufficient to build an effective model.

Something I'm not fully clear about your sampling design is that you seem to keep variants of a given gene in the test set.
Assuming that variants share potentially very similar features, along a common target variable, then my intuition would be to assign each gene to either train or test set.

And for clarification, in a gradient boosting routine, subsample don't refer to splitting data into train and eval/test data. It actually perform a subsample of the train data. The eval data is kept static throughout the training of the collection of trees (it's an optional data that can be provided through the optional deval kwarg. This eval role is typically limited to evaluate best nrounds, and best hyper-parameters. A 3rd data, test data wouldn't be involved at all during the training phase, only kept to assess whether the selected fitted model effectively performs well.

That being said, if there's a desire to have a custom subsampling, my reflex would be to consider a geneID sampling. Ie: at each iteration/tree, extract all variants of a given gene. That could be implemented by selection a sampling of keys from groupby(df, :geneID) assuming a DataFrame based implementation. That would be scenario "2" in my previous comment.

Regarding the implementation, in order to make the information available to subsample, the clean way to procveed I think would involve adding a group ID variable to the cache, similar existing p (prediction) and w (observation weights):

EvoTrees.jl/src/structs.jl

Lines 41 to 42 in c1876dd

	w::Vector{Float32}
	pred::Matrix{Float32}

Cache is set during model instnatiation:

EvoTrees.jl/src/init.jl

Line 106 in c1876dd

cache = CacheBaseCPU{Y,N}(

This could arguably tricky to hack if not already familiar with the codebase though.

[tecosaur]

Is the following a proper depiction of the nature of the dataset: [...]

It is indeed!

Is there some inherent correlation between the number of of variants and the prevalence of the variant class?

Inherently, no. That said, since I pull each class of variants from two different sources, there may be some weak correlation in practice.

If the above condition is respected, my intuition would be that simply splitting the train-eval based on geneID would be sufficient to build an effective model.

Yup, I think so too. The particular scheme I articulated is just what we did to mitigate the complications I listed.

Regarding the implementation, in order to make the information available to subsample, the clean way to proceed I think would involve adding a group ID variable to the cache, similar existing p (prediction) and w (observation weights).

I think that would largely do the trick for us. I suspect grouping + weighting (1/num variants) would be pretty much equivalent to our current sampling scheme, with just one unaddressed factor: the huge class imbalance.

As mentioned, the huge class imbalance is currently handled by only using a (randomly selected) small part of the larger class each bootstrap round, but I'm not sure how to best handle this in a boosted-trees/EvoTrees context.

Assuming that variants share potentially very similar features, along a common target variable, then my intuition would be to assign each gene to either train or test set.

If you squint a little, that's essentially what we're doing here: jut picking a single variant from each gene in the training set, and all variants in the test set.

And for clarification, in a gradient boosting routine, [...]

Thanks! I haven't actually used gradient boosted trees before, and so I still have "traditional" random forest bagging+training rounds in my mind as my base mental model: I've been thinking of the subsampling as an analogue for bootstrapping, but that isn't quite right.

[tecosaur]

Oh @jeremiedb, I wouldn't be surprised if you have just been busy since my last reply, but FWIW I'd be interested to hear more of your thoughts if/when you have the time 🙂

[jeremiedb]

Effectively, I'm a little taken by daily work!

Regarding your problem setup, I find it hard to provide relevant persepctive without having a grasp of what is the really the problem that the model is trying to solve, as I'm missing some understanding of the entites (genes, gene variants, class being predicted).

Regarding the problem setup, I remain confused notably by just picking a single variant from each gene in the training set, and all variants in the test set as I read this as having a some observations (variants) of a given genID into both train and eval datasets.
Here by eval, I'm referring to the deval as found in as a kwarg for fit: https://evovest.github.io/EvoTrees.jl/dev/api/#EvoTrees.fit_evotree

model = EvoTrees.fit(
    config, 
    dtrain;
    target_name,
    feature_names,
    deval)

Can you confirm that tain data is design like:

train_geneIDs = sample(unique(df.geneID))
eval_geneIDs = setdiff(unique(df.geneID), train_geneIDs )
dtrain = df[df.geneID .in train_geneIDs, :]
...

Actually, there would also be a 3rd dataset, dtest, also with its own exclusive set of geneIDs.

Also, am I undersdanding that there are 2 imbalances at play here:

1. Target class imbalance: the model is a logistic regression where the target is Bool where mean(target) is near 0 (or equivalently near 1). By default, I would not perform any fix/hack to adress this situation as the model typically perform just fine even for substantial imbalance. Otherwise, my go-to would be to hack the subsample mechanism happening at each new tree iteration to perform a weighted subsampling that would sample more equal candidates of each class. Then, the model's predictor would need to be scale-back by the sampling bias.
2. Imbalance from the data having a hierarchical structure: individual observations aren't independent, but rather a combinatrion of geneID + variantID. For each geneID, there's a very unequal number of variants. This has 2 implications:
   A. Need to carefully perform the split between dtrain, deval, and dtest (discussed above)
   B. Need to avoid the the over/under importance of geneID specific features caused by the varying number of training observations due to varying number of variants by geneID. I think the most legitimate way to account for this is to compute a weight var like you mentionned that takes the value 1 / num_variants, then specify weight_name = "weight" in fit kwarg.

With the above, you'd then fit like shown in the Titanic tutorial:
https://evovest.github.io/EvoTrees.jl/dev/tutorials/logistic-regression-titanic/

config = EvoTreeRegressor(
  loss=:logloss, 
  nrounds=200, 
  early_stopping_rounds=10,
  eta=0.05, 
  nbins=128, 
  max_depth=5, 
  rowsample=0.5, 
  colsample=0.9)

model = EvoTrees.fit(
    config, dtrain; 
    deval,
    target_name,
    weight_name,
    feature_names,
    print_every_n=10)

[tecosaur]

Thanks for following up Jeremie, I appreciate you lending some of your time to help me work out how I can use EvoTrees with this problem.

I find it hard to provide relevant persepctive without having a grasp of what is the really the problem that the model is trying to solve.

That makes sense, I'm trying to avoid going into an unnecessary level of depth but I'm happy to elaborate where it helps. The problem itself is really just classification with proper handling of the inter-row dependence (within each gene's data) and extreme class imbalance (which reaches peaks around ~10,000 to 1).

I see my description of the current scheme "picking a single variant from each gene in the training set, and all variants in the test set" is a bit ambiguous. I should have said that the test set consists of all variants of genes not in the training set. Your pseudocode looks correct to me in describing the current bootstrap scheme.

The third dataset (dtest) is a surprise to me, but makes sense to me in the context of boosted RF (if I understand it correctly, deval is used in the gradient boosting and dtest is used at the end as an evaluation of the performance independent of the training data).

Also, am I undersdanding that there are 2 imbalances at play here:

Yep, this broadly matches my understanding.

1. One note, when training an RF classifier for this we used stratified sampling to create more balanced training datasets, then calculated performance metrics based on a certain expected class ratio in practice. While I feel this is getting into the weeds, the imbalance in the data we have to train on is different to the imbalance we expect when applying the models.

2. Looks right, nothing to add.

I'll take a look at the Titanic example. From a quick glance look, I'm wondering whether there's a particular reason why you've applied a logistic regression rather than classification? Is that on the expectation that it will behave better in the face of the class imbalance No.1?

[jeremiedb]

For clarificatgion regarding the split of the data into 3 datasets, train, eval and test, it's something not unique to gradient-boosted-trees (like EvoTrees, XGBoost and the likes), but is the canonical aproach common to ML model development that involve hyper-parameter tuning. This is fairly well exposed in: https://en.wikipedia.org/wiki/Training,_validation,_and_test_data_sets

In the setup discussed here, I'd describe it with an added layer of complexity: cross-validations.

1. Split the full dataset into train_full, and test. Here, given the hierarchical nature of the data, sample should be made based on geneID, such all varianats of a geneID belong exclusively to one or the other dataset. This test data won't be involved for the rest of the model training process
2. Split train_full into 2 datasets: train and eval, following the same mutually exclusive geneID allocations. If we are to perform a 5-fold cross-validation, then 5 train / eval pairs are to be built (train-1/eval-1, ..., train-5/eval-5). The each of eval-N an exclusive set of geneIDs, therefore dtrain_full == union(eval-1, ..., eval-5)
3. Hyper-tuning phase: multiple hyper-parameters combination (max_depth, colsample, etc.) are tested. For non-iterative models like random-forest, this can be done by calling fit(learner, dtrain; ...), but for iterative models like EvoTrees, it's a good practice to use an early stopping mechanism by calling fit(learner, train; deval) (this is documented in EvoTrees docs). This is to establish the best number of iterations by training until the evaluation metric stops improving on the eval data. For each of the 5 cross-validation, the metrics are evaluated on the eval-i data associated to the train-i used to fit the model.
4. Once the crvoss-validationh experiments are completed and optimal hyper-parameters identified (including the number of rounds), a final model can be fitted on train-full. This model can then be assessed against the test data.

Note that in the above, all the subsample involved are relevant regardless of the potential subsampling mechanism involved in the ML algorithm itself. It's a valid procedure for GLM, Lasso Regression, etc. In the specific case of EvoTrees, the algorithm performs further subsampling during its training process, so it only affect the specific train-i data called by fit. This is where a particular hacking may be useful in this hierarchical setup, notably by specifying custom weights though weight_name kwargs; and potentially by hacking the subsample used during model training, but I'm not that clear if the later would be particularly usefull. Note that such subsampling could be interpreted similarly to batches used in in neural-network training. Hence, if a custom sampling is deemed usefull in EvoTrees, then it would likely be also be pertinent to build a custom DataLoader to perform hirarchical sampling for a neural network.

[tecosaur]

Thanks for that overview Jeremie, that helps me put together some of what I've read around gradient-boosted RFs. I've also done a bit more reading around, which helps me better contextualise the comments you've made there and previously. I appreciate you taking the time to lay out some of the basics for me :)

It now seems to me that I don't really need much from EvoTrees to be able to apply it to my problem 😀. As long as I do gene-aware splitting into the train/deval/test datasets and weight the data to take into account the varying variant counts and class imbalance I the 4-step process you describe should be viable?

Group-aware subsampling does seem like it could be useful, and from a look around it seems like LightGBM and CatBoost may support this (the docs could be clearer I think): https://lightgbm.readthedocs.io/en/latest/Parameters.html#bagging_by_query, https://catboost.ai/docs/en/concepts/python-reference_pool_set_group_id.

It also looks like using focal loss could be helpful for the sort of extreme imbalance I'm looking at (https://github.com/jhwjhw0123/Imbalance-XGBoost / paper: Imbalance-XGBoost: leveraging weighted and focal losses for binary label-imbalanced classification with XGBoost).

[tecosaur]

@jeremiedb I'm thinking of getting back to this idea some time in the next few months, if you have the opportunity to confirm that my thoughts in my previous comment sound reasonable (or point out anything which doesn't quite seem right), I'd appreciate it 🙂